Changeset 2682


Ignore:
Timestamp:
Apr 10, 2006, 9:39:39 AM (18 years ago)
Author:
nick
Message:

updates to onslow
add print outs to least_squares and mesh

Files:
8 edited

Legend:

Unmodified
Added
Removed

  • inundation/pyvolution/least_squares.py

    r2668 r2682  
    441441        self.mesh = Mesh(vertex_coordinates, triangles,
    442442                         geo_reference = geo)
    443        
     443        print 'mesh.check_integrity'
    444444        self.mesh.check_integrity()
    445 
     445        print 'data_origin'
    446446        self.data_origin = data_origin
    447447
    … …  
    474474                                        data_origin = data_origin,
    475475                                        precrop = precrop)
     476        print 'finished interpolation'
    476477
    477478    def set_point_coordinates(self, point_coordinates,

  • inundation/pyvolution/mesh.py

    r2582 r2682  
    498498
    499499        N = self.number_of_elements
    500 
     500        print 'number of elements:', N
    501501        #Get x,y coordinates for all vertices for all triangles
    502502        V = self.get_vertex_coordinates()
    503 
     503        print 'number of triangle', len(V)
    504504        #Check each triangle
    505505        for i in range(N):
     506           
    506507            x0 = V[i, 0]; y0 = V[i, 1]
    507508            x1 = V[i, 2]; y1 = V[i, 3]
    … …  
    561562        for i in range(N):
    562563            for v in self.triangles[i, :]:
     564                print v
    563565                #Check that all vertices have been registered
    564566                assert self.vertexlist[v] is not None

  • production/onslow_2006/export_results.py

    r2543 r2682  
    11import project, os
     2from os import sep
    23
    34from pyvolution.data_manager import sww2dem
    45from pyvolution.ermapper_grids import read_ermapper_grid
    56
    6 #name = project.outputname
    7 name = project.outputname
     7#name = project.outputname
    88
     9time_dir = "20060328_094539_alpha_0.1"
     10directory = project.outputdir
     11name = directory + time_dir +sep + "source"
     12
     13print'output dir:', name
    914#print 'Which variable do you want to export?'
    1015#which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2, Speed = 3  '  ))
    11 which_var = 2
     16which_var = 4
    1217
    1318if which_var == 0:  # Stage
    … …  
    2934    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-30)'  #Speed
    3035
     36if which_var == 4:  # Elevation
     37    outname = name + '_elevation'
     38    quantityname = 'elevation'  #Elevation
     39
    3140   
    3241sww2dem(name, basename_out = outname,
    3342            quantity = quantityname,
    34             cellsize = 25,       # Trevor would like this at 25
     43            cellsize = 30,       # Trevor would like this at 25
    3544            # define region for viz purposes
    36             easting_min = project.eminviz,
    37             easting_max = project.emaxviz,
    38             northing_min = project.nminviz,
    39             northing_max = project.nmaxviz,       
     45            easting_min = project.e_min_area,
     46            easting_max = project.e_max_area,
     47            northing_min = project.n_min_area,
     48            northing_max = project.n_max_area,       
    4049            reduction = max, #this is because we want max quantityname
    4150            verbose = True,

  • production/onslow_2006/get_timeseries.py

    r2543 r2682  
    77from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    88from pylab import *
    9 from matplotlib.ticker import MultipleLocator, FormatStrFormatter
    10 
    11 swwfile = project.outputname + '.sww'
     9from os import sep
     10# from matplotlib.ticker import MultipleLocator, FormatStrFormatter
     11
     12 
     13
     14time_dir = "20060331_155157"
     15directory = project.outputdir
     16swwfile = directory + time_dir +sep + "source.sww"
     17
     18#swwfile = project.outputname + '.sww'
    1219
    1320def get_gauges_from_file(filename):
    … …  
    2128        fields = line.split(',')
    2229        # my gauge file set up as locationname, easting, northing
    23         location = fields[0]
    24         easting = float(fields[1])
    25         northing = float(fields[2])
     30       
     31        easting = float(fields[0])
     32        northing = float(fields[1])
     33        location = fields[2]
    2634        #z, easting, northing  = redfearn(lat, lon)
    2735        gauges.append([easting, northing])
    … …  
    3139    return gauges, lines, gaugelocation
    3240
    33 gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
     41gauges, lines, gauge_names = get_gauges_from_file(project.gauge_filename)
    3442
    3543print 'number of gauges:   ', len(gauges)
    … …  
    3745#Read model output
    3846quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
     47
    3948f = file_function(swwfile,
    4049                  quantities = quantities,
    … …  
    4756from math import sqrt, atan, degrees
    4857from Numeric import ones
     58
    4959N = len(gauges)
     60
    5061for k, g in enumerate(gauges):
    5162    if k%((N+10)/10)==0: # diagnostics - print 10 lines
    … …  
    7788        uh = f(t, point_id = k)[2]
    7889        vh = f(t, point_id = k)[3]
    79         myloc = locations[k]
     90        myloc = gauge_names[k]
    8091        depth = w-z
    8192
    … …  
    137148           loc='upper right')
    138149    #savefig('Gauge_%d_stage' %k)
    139     savefig('Gauge_%s_stage' %myloc)
    140 
    141     # raw_input('Next')
    142    
     150    #savefig('Gauge_%s_stage' %myloc)
     151    savefig(project.gaugetimeseries+'Gauge_'+myloc+'_stage')
     152
     153    # raw_input('Next')
     154    '''   
    143155    #Momentum plot
    144156    ion()
    … …  
    188200    #savefig('Gauge_%d_speed' %k)
    189201    savefig('Gauge_%s_speed' %myloc)
    190    
     202    '''   
    191203    # raw_input('Next')
    192204

  • production/onslow_2006/project.py

    r2645 r2682  
    4040
    4141#Derive subdirectories and filenames
    42 timedir = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    43 outputdir = home+sep+scenario_dir_name+sep+'output'+sep+timedir+sep
     42time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     43outputtimedir = home+sep+scenario_dir_name+sep+'output'+sep+time+sep
    4444meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
    4545datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
    46 
     46gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep
    4747polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
    4848boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
     49#output dir without time
     50outputdir = home+sep+scenario_dir_name+sep+'output'+sep
     51
     52
    4953print'bound', boundarydir
     54
     55gauge_filename = gaugedir + 'onslow_gauges.xya'
     56
     57gaugetimeseries = gaugedir + 'onslow'
    5058
    5159# boundary source data
    … …  
    6674combined_dem_name = datadir + 'onslow_combined_elevation'
    6775
    68 outputname = outputdir + basename  #Used by post processing
     76outputname = outputtimedir + basename  #Used by post processing
    6977
    7078#!gauge_filename = outputdir + 'onslow_gauges.xya'
    7179#!gauge_outname = outputdir + 'gauges_max_output.xya'
    7280
    73 # clipping region for fine elevation data
    74 eastingmin = 250000
    75 eastingmax = 330000
     81# clipping region to make DEM (pts file) from fine elevation data
     82eastingmin = 240000
     83eastingmax = 340000
    7684northingmin = 7580000
    77 northingmax = 7635000
     85northingmax = 7700000
    7886
    79 south = degminsec2decimal_degrees(-22,00,0)
    80 north = degminsec2decimal_degrees(-21,10,0)
    81 west = degminsec2decimal_degrees(114,30,0)
     87south = degminsec2decimal_degrees(-21,55,0)
     88north = degminsec2decimal_degrees(-20,50,0)
     89west = degminsec2decimal_degrees(114,25,0)
    8290east = degminsec2decimal_degrees(115,30,0)
    83 
     91'''
    8492# region for visualisation
    8593eminviz = 260000
    … …  
    8795nminviz = 7590000
    8896nmaxviz = 7630000
     97'''
     98# region to export
     99
     100e_min_area = 280000
     101e_max_area = 310000
     102n_min_area = 7590000
     103n_max_area = 7630000
    89104
    90105#Georeferencing
    … …  
    93108refzone = 50
    94109
    95 #Main Domain of Onslow: first run NB 21/2/06
    96 d0 = [305000, 7635000]
    97 d1 = [280000, 7635000]
    98 d2 = [250000, 7615000]
    99 d3 = [250000, 7590000]
    100 d4 = [310000, 7580000]
    101 d5 = [330000, 7610000]
    102 d6 = [330000, 7620000]
     110#Updated Main Domain of Onslow: first run NB 6/4/06
     111d0 = [310000, 7690000]
     112d1 = [280000, 7690000]
     113d2 = [270000, 7645000]
     114d3 = [240000, 7625000]
     115d4 = [270000, 7580000]
     116d5 = [300000, 7590000]
     117d6 = [340000, 7610000]
    103118
    104119polyAll = [d0, d1, d2, d3, d4, d5, d6]
    … …  
    122137
    123138poly_thevenard = [j0, j1, j2, j3]
    124 
     139'''
    125140# Direction Is
    126141k0 = [309000, 7619000]
    … …  
    130145
    131146poly_direction = [k0, k1, k2, k3]
    132 
     147'''

  • production/onslow_2006/run_onslow.py

    r2645 r2682  
    33Source data such as elevation and boundary data is assumed to be available in
    44directories specified by project.py
    5 The output sww file is stored in project.outputdir
     5The output sww file is stored in project.outputtimedir
    66
    77The scenario is defined by a triangular mesh created from project.polygon,
    … …  
    4545
    4646# filenames
    47 coarsedemname = project.coarsedemname
     47#coarsedemname = project.coarsedemname
    4848
    4949onshore_dem_name = project.onshore_dem_name
    … …  
    5656
    5757# creates copy of code in output dir if dir doesn't exist
    58 if access(project.outputdir,F_OK) == 0 :
    59     mkdir (project.outputdir)
    60 copy (project.codedirname, project.outputdir + project.codename)
    61 copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py')
     58if access(project.outputtimedir,F_OK) == 0 :
     59    mkdir (project.outputtimedir)
     60copy (project.codedirname, project.outputtimedir + project.codename)
     61copy (project.codedir + 'run_onslow.py', project.outputtimedir + 'run_onslow.py')
    6262print'hello'
     63'''
     64copied_files = False
     65
     66# files to be used
     67files_used = [onshore_dem_name, offshore_points,]
     68
     69if sys.platform != 'win32':   
     70    copied_files = True
     71    for name in file_list:
     72        copy(name, )
     73'''   
    6374#print' most file', project.MOST_dir + project.boundary_basename+'_ha.nc'
    6475#if access(project.MOST_dir + project.boundary_basename+'_ha.nc',F_OK) == 1 :
    6576#    print' most file', project.MOST_dir + project.boundary_basename
    66 
    67 '''
    68 # coarse data
    69 convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
    70 dem2pts(coarsedemname, use_cache=True, verbose=True)
    71 
     77'''
    7278
    7379# fine data (clipping the points file to smaller area)
     80# creates DEM from asc data
    7481convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
     82
     83#creates pts file from DEM
    7584dem2pts(onshore_dem_name,
    7685        easting_min=project.eastingmin,
    … …  
    8089        use_cache=True,
    8190        verbose=True)
    82 '''
    8391
    8492print'create G1'
    … …  
    93101print'export G'
    94102G.new_export_points_file(project.combined_dem_name + '.pts')
    95 
     103'''
    96104
    97105#-------------------------------------------------------------------------------                                 
    … …  
    104112
    105113# original
    106 interior_res = 50000
     114interior_res = 500
    107115interior_regions = [[project.poly_onslow, interior_res],
    108                     [project.poly_thevenard, interior_res],
    109                     [project.poly_direction, interior_res]]
    110                     #[project.testpoly, interior_res]]
     116                    [project.poly_thevenard, interior_res]]
     117
    111118print 'number of interior regions', len(interior_regions)
    112119
    … …  
    115122          project.polyAll,
    116123          {'boundary_tags': {'top': [0], 'topleft': [1],
    117                              'left': [2], 'bottom': [3],
    118                              'bottomright': [4], 'right': [5],
     124                             'topleft1': [2], 'bottomleft': [3],
     125                             'bottom': [4], 'bottomright': [5],
    119126                             'topright':[6]},
    120            'maximum_triangle_area': 1000000,
     127           'maximum_triangle_area': 10000,
    121128           'filename': meshname,           
    122129           'interior_regions': interior_regions},
    … …  
    129136
    130137domain = pmesh_to_domain_instance(meshname, Domain,
    131                                   use_cache = True,
     138                                  use_cache = False,
    132139                                  verbose = True)
    133140
    … …  
    137144
    138145domain.set_name(project.basename)
    139 domain.set_datadir(project.outputdir)
     146domain.set_datadir(project.outputtimedir)
    140147domain.set_quantities_to_be_stored(['stage'])
    141148
    … …  
    190197east = project.east
    191198
     199#note only need to do when an SWW file for the MOST boundary doesn't exist
    192200cache(ferret2sww,
    193201      (source_dir + project.boundary_basename,
    … …  
    228236
    229237domain.set_boundary( {'top': Bf, 'topleft': Bf,
    230                              'left': Br, 'bottom': Br,
    231                              'bottomright': Br, 'right': Br, 'topright': Br} )
    232 
     238                             'topleft1': Bf, 'bottomleft': Bd,
     239                             'bottom': Br, 'bottomright': Br, 'topright': Bd} )
    233240
    234241#-------------------------------------------------------------------------------                                 
    … …  
    238245t0 = time.time()
    239246
    240 for t in domain.evolve(yieldstep = 1000, finaltime = 3000):
     247for t in domain.evolve(yieldstep = 500, finaltime = 3000):
    241248    domain.write_time()
    242249    domain.write_boundary_statistics(tags = 'top')     
    243250
    244 for t in domain.evolve(yieldstep = 100, finaltime = 20000,
    245                        skip_initial_step = True):
     251for t in domain.evolve(yieldstep = 60, finaltime = 15000):
     252#                       ,skip_initial_step = True):
    246253    domain.write_time()
    247254    domain.write_boundary_statistics(tags = 'top')     
Note: See TracChangeset for help on using the changeset viewer.