Changeset 4465 for anuga_work
- Timestamp:
- May 21, 2007, 10:59:13 AM (18 years ago)
- Location:
- anuga_work/production/exmouth_2006
- Files:
-
- 2 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/exmouth_2006/project.py
r4423 r4465 1 # -*- coding: cp1252 -*-2 1 """Common filenames and locations for topographic data, meshes and outputs. 2 Also includes origin for slump scenario. 3 3 """ 4 4 5 from os import sep, environ, getenv, getcwd 6 from os.path import expanduser 5 from os import sep, environ, getenv, getcwd,umask 6 from os.path import expanduser, basename 7 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles 7 8 import sys 9 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees 8 10 from time import localtime, strftime, gmtime 9 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm 11 from anuga.utilities.system_tools import get_user_name, get_host_name 11 12 12 if sys.platform == 'win32': 13 home = getenv('INUNDATIONHOME') 14 user = getenv('USERPROFILE') 13 codename = 'project.py' 15 14 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 20 21 # INUNDATIONHOME is the inundation directory, not the data directory. 22 home += sep +'data' 15 home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir 16 user = get_user_name() 17 host = get_host_name() 18 #needed when running using mpirun, mpirun doesn't inherit umask from .bashrc 19 umask(002) 23 20 24 21 #Making assumptions about the location of scenario data 25 22 state = 'western_australia' 26 scenario_dir_name = 'exmouth_tsunami_scenario_2006' 23 scenario_name = 'exmouth' 24 scenario = 'exmouth_tsunami_scenario' 27 25 28 # onshore data provided by WA DLI 29 onshore_name = 'first' # original 26 #time stuff 27 time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir 28 build_time = time+'_build' 29 run_time = time+'_run' 30 30 31 # 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' 31 #tide = -1.4 32 #tide = 0 33 tide = 1.4 54 34 55 # developed by NM&I 56 coast_name = 'coastline' 35 #Maybe will try to make project a class to allow these parameters to be passed in. 36 alpha = 0.1 37 friction=0.01 38 starttime=3600 39 finaltime=25000 40 setup='trial' 57 41 58 boundary_basename = 'SU-AU' # Mw ? 42 if setup =='trial': 43 print'trial' 44 res_factor=10 45 time_thinning=48 46 yieldstep=240 47 if setup =='basic': 48 print'basic' 49 res_factor=4 50 time_thinning=12 51 yieldstep=120 52 if setup =='final': 53 print'final' 54 res_factor=1 55 time_thinning=1 56 yieldstep=60 59 57 60 #swollen/ all data output 61 basename = 'source' 62 codename = 'project.py' 58 dir_comment='_'+setup+'_'+str(tide)+'_'+str(user) 63 59 64 #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'+sep 67 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep 68 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep 69 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 70 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep 71 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep 72 outputtimedir = outputdir + local_time + sep 73 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 60 onshore_name = 'DLI' 61 onshore_name1 = 'DTED' 62 # offshore 63 offshore_name = 'Exmouth_bathymetry' 64 #offshore_name1 = 'inferred_north' 65 #offshore_name2 = 'inferred_south' 66 coast_name = 'Exmouth_coastline' 74 67 75 gauge_filename = gaugedir + 'exmouth_gauges.csv' 76 community_filename = gaugedir + 'CHINS_v2.csv' 77 community_exmouth = gaugedir + 'community_exmouth.csv' 78 codedir = getcwd()+sep 79 codedirname = codedir + 'project.py' 80 meshname = outputtimedir + 'mesh_' + basename 68 #final topo name 69 combined_name ='exmouth_combined_elevation' 70 #combined_name1 ='exmouth_combined_elevation1' 71 #combined_name_unclipped1 ='exmouth_combined_elevation_unclipped1' 72 combined_small_name = 'exmouth_combined_elevation_small' 81 73 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 74 anuga_dir = home+state+sep+scenario+sep+'anuga'+sep 106 75 107 coast_dem_name = datadir + coast_name 76 topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep 77 topographies_dir = anuga_dir+'topographies'+sep 108 78 109 combined_dem_name = datadir + 'exmouth_combined_elevation' 79 # input topo file location 80 onshore_in_dir_name = topographies_in_dir + onshore_name 81 onshore1_in_dir_name = topographies_in_dir + onshore_name1 82 coast_in_dir_name = topographies_in_dir + coast_name 83 offshore_in_dir_name = topographies_in_dir + offshore_name 84 #offshore_in_dir_name1 = topographies_in_dir + offshore_name1 85 #offshore_in_dir_name2 = topographies_in_dir + offshore_name2 110 86 111 #buildings_filename = gauges_dir + 'exmouth_res_nexis.csv' 87 onshore_dir_name = topographies_dir + onshore_name 88 onshore1_dir_name = topographies_dir + onshore_name1 89 coast_dir_name = topographies_dir + coast_name 90 offshore_dir_name = topographies_dir + offshore_name 91 #offshore_dir_name1 = topographies_dir + offshore_name1 92 #offshore_dir_name2 = topographies_dir + offshore_name2 93 94 #final topo files 95 combined_dir_name = topographies_dir + combined_name 96 #combined_dir_name_unclipped1 = topographies_dir + combined_name_unclipped1 97 #combined_dir_name1 = topographies_dir + combined_name1 98 combined_small_name_dir = topographies_dir + combined_small_name 99 100 meshes_dir = anuga_dir+'meshes'+sep 101 meshes_dir_name = meshes_dir + scenario_name 102 103 polygons_dir = anuga_dir+'polygons'+sep 104 tide_dir = anuga_dir+'tide_data'+sep 105 106 #boundaries locations 107 boundaries_name = 'exmouth_3854_17042007' 108 109 boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'1_10000'+sep 110 boundaries_in_dir_name = boundaries_in_dir + boundaries_name 111 boundaries_dir = anuga_dir+'boundaries'+sep 112 boundaries_dir_name = boundaries_dir + boundaries_name 113 114 #output locations 115 output_dir = anuga_dir+'outputs'+sep 116 output_build_time_dir = output_dir+build_time+sep 117 output_run_time_dir = output_dir +run_time+dir_comment+sep 118 output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing 119 120 #gauges 121 gauge_name = 'exmouth_gauges.csv' 122 gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep 123 gauges_dir_name = gauges_dir + gauge_name 124 125 community_filename = gauges_dir + 'CHINS_v2.csv' 126 community_exmouth = gauges_dir + 'community_exmouth.csv' 127 112 128 buildings_filename_damage_out = 'exmouth_res_nexis_modified.csv' 113 114 community_filename = gaugedir + 'CHINS_v2.csv'115 community_exmouth = gaugedir + 'community_exmouth.csv'116 129 117 130 ############################### 118 131 # Domain definitions 119 132 ############################### 133 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon 120 134 121 # bounding box for clipping MOST/URS output (much bigger than study area) 122 ##south = degminsec2decimal_degrees(-19,0,0) 123 ##north = degminsec2decimal_degrees(-17,15,0) 124 ##west = degminsec2decimal_degrees(120,0,0) 125 ##east = degminsec2decimal_degrees(122,30,0) 126 ## 127 ##d0 = [south, west] 128 ##d1 = [south, east] 129 ##d2 = [north, east] 130 ##d3 = [north, west] 131 ##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3]) 132 ##refzone = zone 135 poly_all = read_polygon(polygons_dir+'extent_zone50a.csv') 133 136 134 # bounding polygon for study area135 poly_all = read_polygon(polygondir+'extent_zone50a.csv')136 137 # plot bounding polygon and make sure BC info surrounds it138 #plot_polygons([poly_all, poly_bc],'boundingpoly',verbose=False)139 137 print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0 140 138 139 res_poly_all = 150000*res_factor 140 141 ############################### 142 # Interior region definitions 143 ############################### 144 ''' 145 poly_0 = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv') 146 res_0 = 20000*res_factor 147 148 poly_1 = read_polygon(polygons_dir+'exmouth_north_coast_inside_extent.csv') 149 res_1 = 5000*res_factor 150 151 poly_2 = read_polygon(polygons_dir+'exmouth_south_coast_inside_extent.csv') 152 res_2 = 5000*res_factor 153 154 poly_3 = read_polygon(polygons_dir+'exmouth_town_pts.csv') 155 res_3 = 2000*res_factor 156 157 poly_4 = read_polygon(polygons_dir+'exmouth_inner_town_pts.csv') 158 res_4 = 500*res_factor 159 160 interior_regions = [[poly_0,res_0],[poly_1,res_1],[poly_2,res_2] 161 ,[poly_3,res_3],[poly_4,res_4]] 162 163 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all) 164 165 print 'min number triangles', trigs_min 166 167 poly_mainland = read_polygon(polygons_dir+'Initial_Condition.csv') 168 ''' 141 169 ################################################################### 142 170 # Clipping regions for export to asc and regions for clipping data … … 148 176 northingmin = 7983427.73 149 177 northingmax = 8032834.52 150 151 ###############################152 # Interior region definitions153 ###############################154 '''155 # exmouth digitized polygons156 poly_exmouth1 = read_polygon(polygondir+'exmouth_Local_Polygon_update.csv')157 poly_exmouth2 = read_polygon(polygondir+'exmouth_Close2_update.csv')158 poly_exmouth3 = read_polygon(polygondir+'exmouth_Coast_update.csv')159 #poly_exmouth4 = read_polygon(polygondir+'Cable_Beach_revised.csv')160 161 plot_polygons([poly_all,poly_exmouth1,poly_exmouth2,poly_exmouth3],'boundingpoly2',verbose=False)162 print 'Area of local polygon', polygon_area(poly_exmouth1)/1000000.0163 print 'Area of close polygon', polygon_area(poly_exmouth2)/1000000.0164 print 'Area of coastal polygon', polygon_area(poly_exmouth3)/1000000.0165 #print 'Area of cable beach polygon', polygon_area(poly_exmouth4)/1000000.0166 167 for i in poly_exmouth3:168 v = is_inside_polygon(i,poly_exmouth1, verbose=False)169 if v == False: print v170 171 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):172 from anuga.utilities.polygon import polygon_area173 174 # TO DO check if any of the regions fall inside one another175 no_triangles = 0.0176 area = polygon_area(bounding_poly)177 for i,j in interior_regions:178 this_area = polygon_area(i)179 no_triangles += this_area/j180 area -= this_area181 print j, this_area/1000000., area/1000000.182 no_triangles += area/remainder_res183 return int(no_triangles/0.7)184 '''
Note: See TracChangeset
for help on using the changeset viewer.