Changeset 4077
- Timestamp:
- Dec 14, 2006, 11:41:51 AM (18 years ago)
- Location:
- anuga_work/production/perth_2006
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/perth_2006/project.py
r3908 r4077 10 10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm 11 11 12 if sys.platform == 'win32': 13 home = getenv('INUNDATIONHOME') 14 user = getenv('USERPROFILE') 12 # file and system info 13 #--------------------------------- 14 codename = 'project.py' 15 15 16 else: 17 home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation') 18 user = getenv('LOGNAME') 19 print 'USER:', user 16 home = getenv('INUNDATIONHOME') #Sandpit's parent dir 17 user = get_user_name() 20 18 21 19 # INUNDATIONHOME is the inundation directory, not the data directory. … … 24 22 #Making assumptions about the location of scenario data 25 23 state = 'western_australia' 24 scenario_name = 'perth' 26 25 scenario_dir_name = 'perth_tsunami_scenario_2006' 27 26 28 27 # onshore data provided by WA DLI 29 onshore_name = '' # original 28 onshore_name = 'perth_dli_ext' # original 29 #island 30 island_name = 'rott_dli_ext' # original 31 island_name1 = 'gard_dli_ext' 32 island_name2 = 'carnac_island_dted' 30 33 31 34 # AHO + DPI data 32 offshore_name1 = 'XY100011610' 33 offshore_name2 = 'XY100011611' 34 offshore_name3 = 'XY100011613' 35 offshore_name4 = 'XY100011614' 36 offshore_name5 = 'XY100011616' 37 offshore_name6 = 'XY100011617' 38 offshore_name7 = 'XY100011618' 39 offshore_name8 = 'XY100011621' 40 offshore_name9 = 'XY100011623' 41 offshore_name10 = 'XY100011745' 42 offshore_name11 = 'XY100011746' 43 offshore_name12 = 'XY100017530' 44 offshore_name13 = 'XY100017532' 45 offshore_name14 = 'XY100017538' 46 offshore_name15 = 'XY100017540' 47 offshore_name16 = 'XYBR66' 48 offshore_name17 = 'XYBR70' 49 offshore_name18 = 'XYBR80' 50 offshore_name19 = 'XYBR88' 51 offshore_name20 = 'XYBR93' 52 offshore_name21 = 'XYBR0110' 53 offshore_name22 = 'XYWADPI' 35 coast_name = 'waterline' 36 offshore_name = 'perth_bathymetry' 54 37 55 # developed by NM&I56 co ast_name = 'coastline'38 #final topo name 39 combined_name ='perth_combined_elevation' 57 40 58 boundary_basename = 'SU-AU' # Mw ? 41 topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep 42 topographies_time_dir = topographies_dir+build_time+sep 59 43 60 #swollen/ all data output 61 basename = 'source' 62 codename = 'project.py' 44 #input topo file location 45 onshore_dir_name = topographies_dir + onshore_name 46 island_dir_name = topographies_dir + island_name 47 island_dir_name1 = topographies_dir + island_name1 48 island_dir_name2 = topographies_dir + island_name2 49 50 coast_dir_name = topographies_dir + coast_name 51 offshore_dir_name = topographies_dir + offshore_name 52 53 #final topo files 54 combined_dir_name = topographies_dir + combined_name 55 combined_time_dir_name = topographies_time_dir + combined_name 56 combined_time_dir_final_name = topographies_time_dir + combined_final_name 57 63 58 64 59 #Derive subdirectories and filenames 65 local_time = strftime('%Y%m%d_%H%M%S',gmtime())66 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep67 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep68 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep69 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep70 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep71 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep72 outputtimedir = outputdir + local_time + sep73 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep74 60 75 gauge_filename = gaugedir + 'perth_gauges.csv' 76 community_filename = gaugedir + 'CHINS_v2.csv' 77 community_perth = gaugedir + 'community_perth.csv' 78 codedir = getcwd()+sep 79 codedirname = codedir + 'project.py' 80 meshname = outputtimedir + 'mesh_' + basename 61 meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep 62 meshes_dir_name = meshes_dir + scenario_name 81 63 82 # Necessary if using point datasets, rather than grid 83 onshore_dem_name = datadir + onshore_name 84 offshore_dem_name1 = datadir + offshore_name1 85 offshore_dem_name2 = datadir + offshore_name2 86 offshore_dem_name3 = datadir + offshore_name3 87 offshore_dem_name4 = datadir + offshore_name4 88 offshore_dem_name5 = datadir + offshore_name5 89 offshore_dem_name6 = datadir + offshore_name6 90 offshore_dem_name7 = datadir + offshore_name7 91 offshore_dem_name8 = datadir + offshore_name8 92 offshore_dem_name9 = datadir + offshore_name9 93 offshore_dem_name10 = datadir + offshore_name10 94 offshore_dem_name11 = datadir + offshore_name11 95 offshore_dem_name12 = datadir + offshore_name12 96 offshore_dem_name13 = datadir + offshore_name13 97 offshore_dem_name14 = datadir + offshore_name14 98 offshore_dem_name15 = datadir + offshore_name15 99 offshore_dem_name16 = datadir + offshore_name16 100 offshore_dem_name17 = datadir + offshore_name17 101 offshore_dem_name18 = datadir + offshore_name18 102 offshore_dem_name19 = datadir + offshore_name19 103 offshore_dem_name20 = datadir + offshore_name20 104 offshore_dem_name21 = datadir + offshore_name21 105 offshore_dem_name22 = datadir + offshore_name22 64 polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep 65 tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep 106 66 107 coast_dem_name = datadir + coast_name 67 #boundaries locations 68 boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep 69 boundaries_in_dir_name = boundaries_in_dir + boundaries_name 70 boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep 71 boundaries_dir_name = boundaries_dir + boundaries_name 72 #boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep 73 #boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing 74 #ideas 75 #boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep 108 76 109 combined_dem_name = datadir + 'perth_combined_elevation' 77 #time stuff 78 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir 79 gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir 80 build_time = time+'_build' 81 run_time = time+'_run' 82 print 'gtime: ', gtime 83 84 #output locations 85 output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep 86 output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep 87 output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep 88 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing 89 90 #gauges 91 gauge_name = 'perth.csv' 92 gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep 93 gauges_dir_name = gauges_dir + gauge_name 94 95 96 97 98 99 100 110 101 111 102 ############################### 112 103 # Domain definitions 113 104 ############################### 105 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 114 106 115 # bounding box for clipping MOST/URS output (much bigger than study area) 116 ##south = degminsec2decimal_degrees(-19,0,0) 117 ##north = degminsec2decimal_degrees(-17,15,0) 118 ##west = degminsec2decimal_degrees(120,0,0) 119 ##east = degminsec2decimal_degrees(122,30,0) 120 ## 121 ##d0 = [south, west] 122 ##d1 = [south, east] 123 ##d2 = [north, east] 124 ##d3 = [north, west] 125 ##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3]) 126 ##refzone = zone 107 bounding_polygon = read_polygon(polygons_dir+'bounding_poly.csv') 127 108 128 # bounding polygon for study area 129 polyAll = read_polygon(polygondir+'extent.csv') 109 refzone = 50 130 110 131 # plot bounding polygon and make sure BC info surrounds it 132 #plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False) 133 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0 111 print "bounding_polygon", bounding_polygon 112 print 'poly area', polygon_area(bounding_polygon)/1000000.0 113 114 115 #Interior regions 116 117 #cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3]) 118 #assert zone == refzone 119 120 poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv') 121 122 poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv') 123 124 poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv') 125 assert zone == refzone 126 134 127 135 128 ################################################################### 136 129 # Clipping regions for export to asc and regions for clipping data 137 130 ################################################################### 138 131 ''' 139 132 # exporting asc grid 140 133 eastingmin = 406215.87 … … 158 151 print 'Area of coastal polygon', polygon_area(poly_perth3)/1000000.0 159 152 #print 'Area of cable beach polygon', polygon_area(poly_perth4)/1000000.0 153 ''' 160 154 161 for i in poly_perth3:162 v = is_inside_polygon(i,poly_perth1, verbose=False)163 if v == False: print v164 165 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):166 from anuga.utilities.polygon import polygon_area167 168 # TO DO check if any of the regions fall inside one another169 no_triangles = 0.0170 area = polygon_area(bounding_poly)171 for i,j in interior_regions:172 this_area = polygon_area(i)173 no_triangles += this_area/j174 area -= this_area175 print j, this_area/1000000., area/1000000.176 no_triangles += area/remainder_res177 return int(no_triangles/0.7)
Note: See TracChangeset
for help on using the changeset viewer.