Changeset 4401


Ignore:
Timestamp:
Apr 18, 2007, 2:44:21 PM (18 years ago)
Author:
nick
Message:

update to dampier

Location:
anuga_work/production/dampier_2006
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/dampier_2006/project.py

    r4373 r4401  
    152152print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
    153153
    154 #res_poly_all = 400000
    155 res_poly_all = 150000
     154res_factor = 1
     155res_poly_all = 150000*res_factor
    156156
    157157###############################
    … …  
    160160
    161161poly_region = read_polygon(polygons_dir+'region.csv')
    162 res_region = 50000
     162res_region = 50000*res_factor
    163163
    164164poly_dampier = read_polygon(polygons_dir+'dampier_town.csv')
    165 #res_dampier = 1000
    166 res_dampier = 500
     165res_dampier = 500*res_factor
    167166
    168167poly_karratha = read_polygon(polygons_dir+'karrathav2.csv')
    169 res_karratha = 15000
     168res_karratha = 15000*res_factor
    170169
    171170poly_karratha_town = read_polygon(polygons_dir+'karratha_townv2.csv')
    172 #res_karratha_town = 1000
    173 res_karratha_town = 500
     171res_karratha_town = 500*res_factor
    174172
    175173poly_facility = read_polygon(polygons_dir+'facility.csv')
    176 res_facility = 1000
     174res_facility = 1000*res_factor
    177175
    178176poly_delambre = read_polygon(polygons_dir+'delambre.csv')
    179 res_delambre = 5000
     177res_delambre = 1000*res_factor
    180178
    181179poly_coast = read_polygon(polygons_dir+'coastpoly.csv')
    182 res_coast = 2500
    183 #res_coast = 10000
     180res_coast = 5000*res_factor
    184181
    185182poly_NWislands = read_polygon(polygons_dir+'nw_islands_area.csv')
    186 res_NWislands = 50000
     183res_NWislands = 50000*res_factor
    187184
    188185poly_island0 = read_polygon(polygons_dir+'island0.csv')
    … …  
    198195res_island0 = res_poly_all
    199196
    200 res_islands = 5000
    201 #res_islands = 15000
     197res_islands = 5000*res_factor
    202198
    203199poly_ref_nw4 = read_polygon(polygons_dir+'ref_nw4.csv')
    … …  
    232228
    233229
    234 ##plot_polygons([poly_dampier,poly_karratha,poly_karratha_town,poly_delambre,
    235 ##                poly_coast,poly_NWislands,poly_island0,poly_island1,poly_island2,
    236 ##                poly_island3,poly_island4,poly_island5,poly_island6,
    237 ##                poly_island7,poly_island8,poly_ref_nw4,poly_ref_nw5,
    238 ##                poly_ref_nw6,poly_ref_nw7,poly_ref_nw8,poly_all],'poly_pic')
     230#plot_polygons([poly_dampier,poly_karratha,poly_karratha_town,poly_delambre,
     231#                poly_coast,poly_NWislands,poly_island0,poly_island1,poly_island2,
     232#                poly_island3,poly_island4,poly_island5,poly_island6,
     233#                poly_island7,poly_island8,poly_ref_nw4,poly_ref_nw5,
     234#                poly_ref_nw6,poly_ref_nw7,poly_ref_nw8,poly_all],'poly_pic')
    239235
    240236interior_regions = [[poly_dampier,res_dampier],
    … …  
    254250                    [poly_karratha,res_karratha],[poly_karratha_town,res_karratha_town]]
    255251
    256 #trigs_min = number_mesh_triangles(interior_regions_test, poly_all, res_poly_all)
    257252trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    258253
    … …  
    263258###################################################################
    264259
    265 poly_bathy = read_polygon(polygons_dir+'mainland_only.csv')
     260poly_mainland = read_polygon(polygons_dir+'mainland_only.csv')
    266261
    267262# exporting asc grid - Dampier
    … …  
    271266n_max_area = 7725000
    272267
    273 # exporting asc grid - Karratha
    274 ##e_min_area =
    275 ##e_max_area =
    276 ##n_min_area =
    277 ##n_max_area =
    278 
    279 """
    280 # used in the CIPMA 2006 scenario
    281 # CIPMA point of interest
    282 cipma_latitude = -20.588456
    283 cipma_longitude = 116.771527
    284 
    285 k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
    286 k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
    287 k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
    288 k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
    289 
    290 cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
    291 assert zone == refzone
    292 
    293 poly_facility = read_polygon(polygons_dir+'facility.csv')
    294 poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
    295 poly_interior = read_polygon(polygons_dir+'interior.csv')
    296 poly_coast = read_polygon(polygons_dir+'coast_final.csv')
    297 clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
    298 clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
    299 clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
    300 clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
    301 
    302 # exporting asc grid
    303 e_min_area = 474000
    304 e_max_area = 480000
    305 n_min_area = 7719000
    306 n_max_area = 7725000
    307 
    308 # used in the original 2005 scenario
    309 #Interior regions
    310 karratha_south = degminsec2decimal_degrees(-20,44,0)
    311 karratha_north = degminsec2decimal_degrees(-20,42,0)
    312 karratha_west = degminsec2decimal_degrees(116,48,0)
    313 karratha_east = degminsec2decimal_degrees(116,53,30)
    314 
    315 k0 = [karratha_south, karratha_west]
    316 k1 = [karratha_south, karratha_east]
    317 k2 = [karratha_north, karratha_east]
    318 k3 = [karratha_north, karratha_west]   
    319 
    320 karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
    321 assert zone == refzone
    322 
    323 #Interior regions
    324 dampier_south = degminsec2decimal_degrees(-20,40,0)
    325 dampier_north = degminsec2decimal_degrees(-20,38,10)
    326 dampier_west = degminsec2decimal_degrees(116,43,0)
    327 dampier_east = degminsec2decimal_degrees(116,45,0)
    328 
    329 d0 = [dampier_south, dampier_west]
    330 d1 = [dampier_south, dampier_east]
    331 d2 = [dampier_north, dampier_east]
    332 d3 = [dampier_north, dampier_west]   
    333 
    334 dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
    335 assert zone == refzone
    336 
    337 #Interior regions
    338 refinery_south = degminsec2decimal_degrees(-20,37,50)
    339 refinery_north = degminsec2decimal_degrees(-20,36,0)
    340 refinery_west = degminsec2decimal_degrees(116,44,0)
    341 refinery_east = degminsec2decimal_degrees(116,46,10)
    342 
    343 d0 = [refinery_south, refinery_west]
    344 d1 = [refinery_south, refinery_east]
    345 d2 = [refinery_north, refinery_east]
    346 d3 = [refinery_north, refinery_west]   
    347 
    348 refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
    349 assert zone == refzone
    350 
    351 #Interior region around 468899, 7715177:
    352 #lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
    353 
    354 point_south = degminsec2decimal_degrees(-20,39,46)
    355 point_north = degminsec2decimal_degrees(-20,39,42)
    356 point_west = degminsec2decimal_degrees(116,42,0)
    357 point_east = degminsec2decimal_degrees(116,42,10)
    358 
    359 d0 = [point_south, point_west]
    360 d1 = [point_south, point_east]
    361 d2 = [point_north, point_east]
    362 d3 = [point_north, point_west]   
    363 
    364 point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
    365 assert zone == refzone
    366 
    367 #Neils areas around interesting points
    368 neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
    369                 degminsec2decimal_degrees(116,45,18)]
    370 neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
    371                 degminsec2decimal_degrees(116,46,18)]
    372 neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
    373                 degminsec2decimal_degrees(116,47,17)]
    374 neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
    375                 degminsec2decimal_degrees(116,46,17)]
    376 
    377 neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
    378                                                   neil1_point2,
    379                                                   neil1_point3,
    380                                                   neil1_point4])
    381 assert zone == refzone
    382 
    383 neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
    384                 degminsec2decimal_degrees(116,41,33)]
    385 neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
    386                 degminsec2decimal_degrees(116,42,13)]
    387 neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
    388                 degminsec2decimal_degrees(116,43,49)]
    389 neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
    390                 degminsec2decimal_degrees(116,43,9)]
    391 
    392 neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
    393                                                   neil2_point2,
    394                                                   neil2_point3,
    395                                                   neil2_point4])
    396 assert zone == refzone
    397 
    398 #Withnell bay
    399 wb_point1 = [degminsec2decimal_degrees(-20,35,34),
    400                 degminsec2decimal_degrees(116,45,18)]
    401 wb_point2 = [degminsec2decimal_degrees(-20,36,15),
    402                 degminsec2decimal_degrees(116,46,18)]
    403 wb_point3 = [degminsec2decimal_degrees(-20,35,9),
    404                 degminsec2decimal_degrees(116,47,17)]
    405 wb_point4 = [degminsec2decimal_degrees(-20,34,26),
    406                 degminsec2decimal_degrees(116,46,17)]
    407 
    408 wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
    409                                                wb_point3, wb_point4])
    410 assert zone == refzone
    411 
    412 #Larger Withnell bay
    413 lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
    414                 degminsec2decimal_degrees(116,42,00)]
    415 lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
    416                 degminsec2decimal_degrees(116,46,50)]
    417 lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
    418                 degminsec2decimal_degrees(116,47,39)]
    419 lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
    420                 degminsec2decimal_degrees(116,42,50)]
    421 
    422 lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
    423                                                 lwb_point3, lwb_point4])
    424                                                      
    425 assert zone == refzone
    426 """

  • anuga_work/production/dampier_2006/run_boundary_points.py

    r4364 r4401  
    3030print "poly", poly
    3131print "project.boundaries_dir_name5",project.boundaries_dir_name5
    32 URS_points_needed_to_file(project.boundaries_dir_name5,poly, export_csv=True, zone=50)
     32URS_points_needed_to_file('i',poly, export_csv=True, zone=50)
     33#URS_points_needed_to_file(project.boundaries_dir_name5,poly, export_csv=True, zone=50)
    3334#URS_points_needed_to_file(project.boundaries_dir_name5,poly, export_csv=True)
    3435#urs_ungridded2sww(project.boundaries_in_dir_name2, project.boundaries_dir_name5, mean_stage= project.tide,  verbose=True)

  • anuga_work/production/dampier_2006/run_dampier.py

    r4373 r4401  
    7575# resolutions (maximal area of per triangle) for each polygon
    7676#--------------------------------------------------------------------------
    77 '''
    78 poly = [[0,0],[0,100],[100,100],[100,0]]
    79 
    80 create_mesh_from_regions(poly,
    81                              boundary_tags={'back': [0], 'side': [1,3],
    82                                             'ocean': [2]},
    83                              maximum_triangle_area=1,
    84                              interior_regions=None,
    85                              filename=meshes_dir_name,
    86                              use_cache=True,
    87                              verbose=True)
    88 
    89 sys.exit()
    90 '''
     77
    9178if myid == 0:
    9279   
    9380    print 'start create mesh from regions'
    94 #    print 'project.res_poly_all',project.res_poly_all
    95 #    print 'project.interior_regions_test',project.interior_regions_test
    96 #    print 'meshes_dir_name',meshes_dir_name
     81
    9782    create_mesh_from_regions(project.poly_all,
    9883                             boundary_tags={'back': [2,3], 'side': [0, 1, 4],
    … …  
    10186                             interior_regions=project.interior_regions,
    10287                             filename=meshes_dir_name,
    103                              use_cache=True,
    104                              verbose=True)
    105 
    106 # to sync all processors are ready
    107 barrier()
    108 '''
    109 create_mesh_from_regions(project.poly_all,
    110                              boundary_tags={'back': [2,3], 'side': [0, 1, 4],
    111                                             'ocean': [5]},
    112                              maximum_triangle_area=project.res_poly_all,
    113 #                             interior_regions=project.interior_regions,
    114                              interior_regions=None,
    115                              filename=meshes_dir_name,
    11688                             use_cache=False,
    11789                             verbose=True)
    118 '''
     90
     91# to sync all processors are ready
     92barrier()
     93
    11994#-------------------------------------------------------------------------
    12095# Setup computational domain
    … …  
    136111
    137112from anuga.shallow_water.data_manager import urs2sww
    138 '''
    139 # put above distribute
    140 print 'boundary file is: ',project.boundaries_dir_name
    141 from caching import cache
    142 if myid == 0:
    143     cache(urs2sww,
    144           (project.boundaries_in_dir_name,
    145            project.boundaries_dir_name),
    146           {'verbose': True,
    147            'minlat': project.south_boundary,
    148            'maxlat': project.north_boundary,
    149            'minlon': project.west_boundary,
    150            'maxlon': project.east_boundary,
    151            'mint': 0, 'maxt': 35100,
    152            'origin': domain.geo_reference.get_origin(),
    153            'mean_stage': project.tide,
    154 #           'zscale': 1,                 #Enhance tsunami
    155            'fail_on_NaN': False},
    156            verbose = True,
    157            )
    158 barrier()
    159 '''
     113
    160114
    161115#-------------------------------------------------------------------------
    … …  
    168122    from polygon import *
    169123    #following sets the stage/water to be offcoast only
    170     IC = Polygon_function( [(project.poly_bathy, 0.)], default = tide)
     124    IC = Polygon_function( [(project.poly_mainland, 0.)], default = tide)
     125#    IC = Polygon_function( [(project.poly_mainland, 0.)], default = 200)
    171126    domain.set_quantity('stage', IC)
    172127    domain.set_quantity('friction', 0.01)
    … …  
    177132# MUST USE TXT FILES FOR CACHING TO WORK!
    178133                    filename = project.combined_dir_name + '.txt',
     134#                    filename = project.combined_smallest_dir_name + '.txt',
    179135                    use_cache = True,
    180136                    verbose = True,
    … …  
    211167print 'Available boundary tags', domain.get_boundary_tags()
    212168print 'domain id', id(domain)
    213 print 'Reading Boundary file',project.boundaries_dir_name4 + '.sww'
    214 #Bf = File_boundary(boundaries_dir_name + '.sww',
    215 #Bf = File_boundary(project.boundaries_dir_name3 + '.sww',
    216 #                  domain, time_thinning=24, use_cache=True, verbose=True)
     169#print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
     170
    217171Bf = Field_boundary(project.boundaries_dir_namea + '.sww',
    218                   domain, time_thinning=12, mean_stage=tide, use_cache=True, verbose=True)
     172                    domain, time_thinning=12, mean_stage=tide,
     173                    use_cache=True, verbose=True)
    219174
    220175print 'finished reading boundary file'
    … …  
    238193t0 = time.time()
    239194
    240 for t in domain.evolve(yieldstep = 60, finaltime = 29990):
     195for t in domain.evolve(yieldstep = 60, finaltime = 29950):
    241196    domain.write_time()
     197#    domain.write_time(track_speeds=True)
    242198    domain.write_boundary_statistics(tags = 'ocean')
    243199   
    244     domain.write_time(track_speeds=True)
     200#    domain.write_time(track_speeds=True)
    245201
    246202#for t in domain.evolve(yieldstep = 120, finaltime = 9000):
Note: See TracChangeset for help on using the changeset viewer.