Changeset 3881
- Timestamp:
- Oct 30, 2006, 9:21:02 AM (18 years ago)
- Location:
- anuga_work/production/broome_2006
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/broome_2006/project.py
r3788 r3881 8 8 from time import localtime, strftime, gmtime 9 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_utm10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm 11 11 12 12 if sys.platform == 'win32': … … 27 27 28 28 # onshore data provided by WA DLI 29 onshore_name = ' ' # original29 onshore_name = 'dted2_z51' # original 30 30 31 # offshore data provided by WA DPI 32 offshore_name_dpi1 = '' 33 34 # AHO data 35 offshore_name1 = 'xy100003760' 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' 36 54 37 55 # developed by NM&I 38 coast_name = 'coastline _points'56 coast_name = 'coastline' 39 57 40 58 boundary_basename = 'SU-AU' # Mw ? … … 63 81 # Necessary if using point datasets, rather than grid 64 82 onshore_dem_name = datadir + onshore_name 65 offshore_dem_name_local1 = datadir + offshore_name_dpi1 66 offshore_dem_name_aho1 = datadir + offshore_name1 83 offshore_dem_name1 = datadir + offshore_name1 84 offshore_dem_name2 = datadir + offshore_name2 85 offshore_dem_name3 = datadir + offshore_name3 86 offshore_dem_name4 = datadir + offshore_name4 87 offshore_dem_name5 = datadir + offshore_name5 88 offshore_dem_name6 = datadir + offshore_name6 89 offshore_dem_name7 = datadir + offshore_name7 90 offshore_dem_name8 = datadir + offshore_name8 91 offshore_dem_name9 = datadir + offshore_name9 92 offshore_dem_name10 = datadir + offshore_name10 93 offshore_dem_name11 = datadir + offshore_name11 94 offshore_dem_name12 = datadir + offshore_name12 95 offshore_dem_name13 = datadir + offshore_name13 96 offshore_dem_name14 = datadir + offshore_name14 97 offshore_dem_name15 = datadir + offshore_name15 98 offshore_dem_name16 = datadir + offshore_name16 99 offshore_dem_name17 = datadir + offshore_name17 100 offshore_dem_name18 = datadir + offshore_name18 101 offshore_dem_name19 = datadir + offshore_name19 102 offshore_dem_name20 = datadir + offshore_name20 103 offshore_dem_name21 = datadir + offshore_name21 104 offshore_dem_name22 = datadir + offshore_name22 67 105 68 106 coast_dem_name = datadir + coast_name … … 75 113 76 114 # bounding box for clipping MOST/URS output (much bigger than study area) 77 south = degminsec2decimal_degrees(-19,0,0)78 north = degminsec2decimal_degrees(-17,15,0)79 west = degminsec2decimal_degrees(120,0,0)80 east = degminsec2decimal_degrees(122,30,0)81 82 d0 = [south, west]83 d1 = [south, east]84 d2 = [north, east]85 d3 = [north, west]86 poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])87 refzone = zone115 ##south = degminsec2decimal_degrees(-19,0,0) 116 ##north = degminsec2decimal_degrees(-17,15,0) 117 ##west = degminsec2decimal_degrees(120,0,0) 118 ##east = degminsec2decimal_degrees(122,30,0) 119 ## 120 ##d0 = [south, west] 121 ##d1 = [south, east] 122 ##d2 = [north, east] 123 ##d3 = [north, west] 124 ##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3]) 125 ##refzone = zone 88 126 89 127 # bounding polygon for study area … … 91 129 92 130 # plot bounding polygon and make sure BC info surrounds it 93 plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)131 #plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False) 94 132 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0 95 133 … … 98 136 ################################################################### 99 137 100 # clipping 20m onshore data set 101 #eastingmin = 520000 102 #eastingmax = 536000 103 #northingmin = 5245000 104 #northingmax = 5260000 138 # exporting asc grid 139 eastingmin = 406215.87 140 eastingmax = 440208.78 141 northingmin = 7983427.73 142 northingmax = 8032834.52 143 144 145 146 105 147 106 148 ############################### … … 112 154 poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv') 113 155 poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv') 114 poly_broome4 = read_polygon(polygondir+'Cable_Beach .csv')156 poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv') 115 157 116 158 plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False) … … 123 165 v = is_inside_polygon(i,poly_broome1, verbose=False) 124 166 if v == False: print v 167 168 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res): 169 from anuga.utilities.polygon import polygon_area 170 171 # TO DO check if any of the regions fall inside one another 172 no_triangles = 0.0 173 area = polygon_area(bounding_poly) 174 for i,j in interior_regions: 175 this_area = polygon_area(i) 176 no_triangles += this_area/j 177 area -= this_area 178 print j, this_area/1000000., area/1000000. 179 no_triangles += area/remainder_res 180 return int(no_triangles/0.7) -
anuga_work/production/broome_2006/run_broome.py
r3788 r3881 10 10 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 11 11 """ 12 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res): 13 from anuga.utilities.polygon import polygon_area 14 15 # check if any of the regions fall inside one another 16 no_triangles = 0.0 17 area = polygon_area(bounding_poly) 18 for i,j in interior_regions: 19 this_area = polygon_area(i) 20 no_triangles += this_area/j 21 area -= this_area 22 print j, this_area/1000000., area/1000000. 23 no_triangles += area/remainder_res 24 return int(no_triangles/0.7) 12 25 13 #------------------------------------------------------------------------------- 26 14 # Import necessary modules … … 74 62 75 63 # filenames 76 #onshore_dem_name = project.onshore_dem_name 77 #coast_points = project.coast_dem_name 78 #offshore_points = project.offshore_dem_name 64 onshore_dem_name = project.onshore_dem_name 65 coast_points = project.coast_dem_name 79 66 meshname = project.meshname+'.msh' 80 #source_dir = project.boundarydir 81 82 copied_files = False 83 ''' 67 84 68 # creates DEM from asc data 85 69 convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True) 86 70 87 71 #creates pts file for onshore DEM 88 dem2pts(onshore_dem_name, 89 easting_min=project.eastingmin, 90 easting_max=project.eastingmax, 91 northing_min=project.northingmin, 92 northing_max= project.northingmax, 93 use_cache=True, 94 verbose=True) 72 dem2pts(onshore_dem_name, use_cache=True, verbose=True) 95 73 96 74 print 'create offshore' 97 G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya') 75 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\ 76 Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\ 77 Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\ 78 Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\ 79 Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\ 80 Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\ 81 Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\ 82 Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\ 83 Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\ 84 Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\ 85 Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\ 86 Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\ 87 Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\ 88 Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\ 89 Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\ 90 Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\ 91 Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\ 92 Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\ 93 Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\ 94 Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\ 95 Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\ 96 Geospatial_data(file_name = project.offshore_dem_name22 + '.xya') 98 97 print 'create onshore' 99 98 G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') … … 101 100 G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya') 102 101 print 'add' 103 G = G1 + G2 + G3 + G4102 G = G1 + G2 + G3 104 103 print 'export points' 105 104 G.export_points_file(project.combined_dem_name + '.pts') 106 ''' 105 107 106 #---------------------------------------------------------------------------- 108 107 # Create the triangular mesh based on overall clipping polygon with a tagged … … 112 111 113 112 from anuga.pmesh.mesh_interface import create_mesh_from_regions 114 remainder_res = 100000115 local_res = 15000116 broome_res = 2500113 remainder_res = 750000 114 local_res = 25000 115 broome_res = 5000 117 116 coast_res = 500 118 117 beach_res = 500 … … 122 121 [project.poly_broome4, beach_res]] 123 122 124 print 'number of interior regions', len(interior_regions) 123 from project import number_mesh_triangles 125 124 print 'estimate of number of triangles', \ 126 125 number_mesh_triangles(interior_regions, project.polyAll, remainder_res) 127 128 import sys; sys.exit()129 126 130 127 from caching import cache … … 139 136 verbose = True, evaluate=False) 140 137 print 'created mesh' 141 import sys; sys.exit() 138 142 139 #------------------------------------------------------------------------------- 143 140 # Setup computational domain … … 161 158 domain.set_quantity('stage', tide) 162 159 domain.set_quantity('friction', 0.0) 163 domain.set_quantity('elevation', 0.0,164 #filename = project.combined_dem_name + '.pts',160 domain.set_quantity('elevation', 161 filename = project.combined_dem_name + '.pts', 165 162 use_cache = True, 166 163 verbose = True, … … 221 218 t0 = time.time() 222 219 223 for t in domain.evolve(yieldstep = 240, finaltime = 6800):220 for t in domain.evolve(yieldstep = 240, finaltime = 480): 224 221 domain.write_time() 225 222 domain.write_boundary_statistics(tags = 'e14') 226 227 for t in domain.evolve(yieldstep = 30, finaltime = 9000228 ,skip_initial_step = True):229 domain.write_time()230 domain.write_boundary_statistics(tags = 'e14')231 232 for t in domain.evolve(yieldstep = 240, finaltime = 15000233 ,skip_initial_step = True):234 domain.write_time()235 domain.write_boundary_statistics(tags = 'e14')236 223 237 224 print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.