Changeset 3788
- Timestamp:
- Oct 16, 2006, 2:12:04 PM (18 years ago)
- Location:
- anuga_work/production
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/broome_2006/project.py
r3769 r3788 1 # -*- coding: cp1252 -*- 1 2 """Common filenames and locations for topographic data, meshes and outputs. 2 3 """ 3 4 4 5 from os import sep, environ, getenv, getcwd 5 from os.path import expanduser, basename 6 #from anuga.utilities.polygon import read_polygon 6 from os.path import expanduser 7 7 import sys 8 from anuga.pmesh.create_mesh import convert_from_latlon_to_utm 9 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees 10 from time import localtime, strftime 11 from anuga.geospatial_data.geospatial_data import * 8 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 12 if sys.platform == 'win32': 13 home = getenv('INUNDATIONHOME') 14 user = getenv('USERPROFILE') 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 20 21 # INUNDATIONHOME is the inundation directory, not the data directory. 22 home += sep +'data' 12 23 13 24 #Making assumptions about the location of scenario data … … 15 26 scenario_dir_name = 'broome_tsunami_scenario_2006' 16 27 17 # all data to be delivered by National Mapping 18 # onshore data from 30m DTED level 2 19 onshore_name_dted = 'broome_onshore_30m_dted' 20 onshore_name_dli = 'broome_onshore_20m_dli' 28 # onshore data provided by WA DLI 29 onshore_name = '' # original 21 30 22 # offshore data from GA digitised charts23 offshore_name 1 = 'broome_offshore_points'31 # offshore data provided by WA DPI 32 offshore_name_dpi1 = '' 24 33 25 # offshore data from AHO fairsheets26 offshore_name 2 = 'broome_offshore_points_fairsheet'34 # AHO data 35 offshore_name1 = 'xy100003760' 27 36 28 # coastline developed from aerial photography and 1.5m DLI contour29 coast_name = ' broome_coastline_points'37 # developed by NM&I 38 coast_name = 'coastline_points' 30 39 31 boundary_basename = 'SU-AU _clip'40 boundary_basename = 'SU-AU' # Mw ? 32 41 33 42 #swollen/ all data output … … 35 44 codename = 'project.py' 36 45 37 if sys.platform == 'win32':38 home = getenv('INUNDATIONHOME') #Sandpit's parent dir39 user = getenv('USERPROFILE')40 else:41 # update to perlite 242 home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'2'+sep+'cit'+sep+'inundation'+sep+'data')43 user = getenv('LOGNAME')44 45 # INUNDATIONHOME is the inundation directory, not the data directory.46 home += sep +'data'47 48 46 #Derive subdirectories and filenames 49 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir 50 outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep 51 47 local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 52 48 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep 53 49 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep 54 50 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep 55 51 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 56 boundarydir = home+sep+s tate+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep52 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep 57 53 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep 58 tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep 54 outputtimedir = outputdir + local_time + sep 55 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep 59 56 60 gauge_filename = gaugedir + 'gauge_location_broome.csv' 61 buildings_filename = gaugedir + 'broome_res.csv' 62 buildings_filename_out = 'broome_res_modified.csv' 63 buildings_filename_damage_out = 'broome_res_modified_damage.csv' 64 community_filename = gaugedir + 'CHINS_v2.csv' 65 community_scenario = gaugedir + 'community_pt_hedland.csv' 66 tidal_filename = tidedir + 'pt_hedland_tide.txt' 57 gauge_filename = gaugedir + 'broome_gauges.csv' 67 58 68 meshname = meshdir + basename 69 onshore_dem_name = datadir + onshore_name_dli 70 offshore_dem_name1 = datadir + offshore_name1 71 offshore_dem_name2 = datadir + offshore_name2 59 codedir = getcwd()+sep 60 codedirname = codedir + 'project.py' 61 meshname = outputtimedir + 'mesh_' + basename 62 63 # Necessary if using point datasets, rather than grid 64 onshore_dem_name = datadir + onshore_name 65 offshore_dem_name_local1 = datadir + offshore_name_dpi1 66 offshore_dem_name_aho1 = datadir + offshore_name1 67 72 68 coast_dem_name = datadir + coast_name 73 combined_dem_name = datadir + 'broome_combined_elevation'74 outputname = outputtimedir + basename #Used by post processing75 69 76 # for ferret2sww 77 south = degminsec2decimal_degrees(-20,30,0) 78 north = degminsec2decimal_degrees(-17,10,0) 79 west = degminsec2decimal_degrees(117,00,0) 80 east = degminsec2decimal_degrees(120,00,0) 70 combined_dem_name = datadir + 'broome_combined_elevation' 81 71 82 # region to export (used from export_results.py) 83 e_min_area = 648000 84 e_max_area = 675000 85 n_min_area = 7745000 86 n_max_area = 7761000 72 ############################### 73 # Domain definitions 74 ############################### 87 75 88 export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],[e_max_area,n_max_area],[e_max_area,n_min_area]] 89 90 refzone = 50 76 # 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) 91 81 92 from anuga.coordinate_transforms.redfearn import redfearn 93 # boundary up to 50 m contour 94 lat1_50 = degminsec2decimal_degrees(-19,20,0) 95 lat2_50 = degminsec2decimal_degrees(-19,30,0) 96 lat3_50 = degminsec2decimal_degrees(-19,45,0) 97 lon1_50 = degminsec2decimal_degrees(119,05,0) 98 lon2_50 = degminsec2decimal_degrees(118,20,0) 99 lon3_50 = degminsec2decimal_degrees(117,45,0) 100 z, easting, northing = redfearn(lat1_50, lon1_50) 101 d0_50 = [easting, northing] 102 z, easting, northing = redfearn(lat2_50, lon2_50) 103 d1_50 = [easting, northing] 104 z, easting, northing= redfearn(lat3_50, lon3_50) 105 d2_50 = [easting, northing] 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 = zone 106 88 107 d4_50 = [285000, 7585000] 108 d6_50 = [330000, 7605000] 109 #bounding_poly50 = [p0_50, p1_50, p2_50, d6_50, d5, d4_50] 89 # bounding polygon for study area 90 polyAll = read_polygon(polygondir+'extent.csv') 110 91 111 d0 = [763852.0, 7934358.0] 112 d1 = [710987.0, 7925797.0] 113 d2 = [658264.0, 7926314.0] 114 d3 = [552686.0, 7871580.0] 115 #d4 = [604415.81, 7733013.56] 116 d4 = [638000.0, 7733013.56] 117 #d5 = [656561.15, 7732615.11] 118 d5 = [662000.0, 7732615.11] 119 #d6 = [708940.32, 7750510.33] 120 d6 = [690000.0, 7740510.33] 121 #polyAll = [d0, d1, d2, d3, d4, d5, d6] 122 #polyAll = [d0_50, d1_50, d2_50, d4, d5, d6] 123 # from Hamish 124 h0=[629262.17, 7747205.47] 125 h1=[552686.00, 7871579.99] #d3 126 h2=[658264.00, 7926314.00] #d2 127 h3=[710986.99, 7925796.99] #d1 128 h4=[763851.99, 7934357.99] #d0 129 h5=[701485.21, 7770656.86] 130 h6=[698273.75, 7762227.38] 131 h7=[698194.23, 7762018.65] 132 h8=[691627.41, 7744781.98] 133 h9=[679220.75, 7743604.59] 134 h10=[653512.59, 7740528.56] 135 h11=[634777.71, 7738247.17] 136 h12=[629443.86, 7746910.37] 137 h13=[629396.84, 7746986.75] 138 h14=[629352.32, 7747059.06] 139 h15=[629276.24, 7747182.63] 140 h16=[629262.17, 7747205.47] #repeat of h0 141 # using Hamish's new bounding polygon 142 #polyAll = [d0_50, d1_50, d2_50, h16,h15,h14,h13,h12,h11,h10,h9,h8,h7,h6,h5] 143 polyAll = [d0_50, d1_50, d2_50, h16,h11,h8,h6, h5] 92 # plot bounding polygon and make sure BC info surrounds it 93 plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False) 94 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0 144 95 145 #Interior region - Broome centre 146 i0 = [668000, 7757000] 147 i1 = [659000, 7755000] 148 i2 = [660000, 7749000] 149 i3 = [667000, 7746000] 150 i4 = [678000, 7751000] 96 ################################################################### 97 # Clipping regions for export to asc and regions for clipping data 98 ################################################################### 151 99 152 poly_broome = [i0, i1, i2, i3, i4] 100 # clipping 20m onshore data set 101 #eastingmin = 520000 102 #eastingmax = 536000 103 #northingmin = 5245000 104 #northingmax = 5260000 153 105 154 #Are there other significant features? 155 j0 = [670000, 7760000] 156 j1 = [633000, 7745000] 157 j2 = [665000, 7743000] 158 j3 = [690000, 7755000] 106 ############################### 107 # Interior region definitions 108 ############################### 159 109 160 poly_region = [j0, j1, j2, j3] 110 # broome digitized polygons 111 poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon.csv') 112 poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv') 113 poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv') 114 poly_broome4 = read_polygon(polygondir+'Cable_Beach.csv') 115 116 plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False) 117 print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0 118 print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0 119 print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0 120 print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0 121 122 for i in poly_broome3: 123 v = is_inside_polygon(i,poly_broome1, verbose=False) 124 if v == False: print v -
anuga_work/production/broome_2006/run_broome.py
r3535 r3788 6 6 7 7 The scenario is defined by a triangular mesh created from project.polygon, 8 the elevation data and a simulated submarine landslide.8 the elevation data and a tsunami wave generated by MOST. 9 9 10 10 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006 11 11 """ 12 #-------------------------------------------------------------------------------# Import necessary modules 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) 25 #------------------------------------------------------------------------------- 26 # Import necessary modules 13 27 #------------------------------------------------------------------------------- 14 28 15 29 # Standard modules 16 from os import sep 17 from os.path import dirname, basename 30 import os 18 31 import time 19 20 # Related major packages21 from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \22 Dirichlet_boundary, Time_boundary, File_boundary23 from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts24 from anuga.pyvolution.combine_pts import combine_rectangular_points_files25 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance26 32 from shutil import copy 27 33 from os import mkdir, access, F_OK 34 import sys 35 36 # Related major packages 37 from anuga.shallow_water import Domain, Reflective_boundary, \ 38 Dirichlet_boundary, Time_boundary, File_boundary 39 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts 40 from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 28 41 from anuga.geospatial_data.geospatial_data import * 29 import sys 30 from anuga.pyvolution.util import Screen_Catcher 42 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher 31 43 32 44 # Application specific imports … … 41 53 if access(project.outputtimedir,F_OK) == 0 : 42 54 mkdir (project.outputtimedir) 43 copy ( dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')44 copy ( __file__, project.outputtimedir + basename(__file__))45 print 'project.outputtimedir',project.outputtimedir46 47 # 55 copy (project.codedirname, project.outputtimedir + project.codename) 56 copy (project.codedir + 'run_broome.py', project.outputtimedir + 'run_broome.py') 57 print'output dir', project.outputtimedir 58 59 #normal screen output is stored in 48 60 screen_output_name = project.outputtimedir + "screen_output.txt" 49 61 screen_error_name = project.outputtimedir + "screen_error.txt" 50 62 51 # 63 #used to catch screen output to file 52 64 sys.stdout = Screen_Catcher(screen_output_name) 53 65 sys.stderr = Screen_Catcher(screen_error_name) 66 54 67 print 'USER: ', project.user 55 68 … … 58 71 # 59 72 # Convert ASC 2 DEM 2 PTS using source data and store result in source data 60 # Do for coarse and fine data61 # Fine pts file to be clipped to area of interest62 73 #------------------------------------------------------------------------------- 63 74 64 75 # filenames 76 #onshore_dem_name = project.onshore_dem_name 77 #coast_points = project.coast_dem_name 78 #offshore_points = project.offshore_dem_name 65 79 meshname = project.meshname+'.msh' 66 source_dir = project.boundarydir 67 68 # creates DEM from asc data 69 convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True) 70 71 #creates pts file from DEM 72 dem2pts(project.onshore_dem_name, 80 #source_dir = project.boundarydir 81 82 copied_files = False 83 ''' 84 # creates DEM from asc data 85 convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True) 86 87 #creates pts file for onshore DEM 88 dem2pts(onshore_dem_name, 73 89 easting_min=project.eastingmin, 74 90 easting_max=project.eastingmax, 75 91 northing_min=project.northingmin, 76 92 northing_max= project.northingmax, 77 use_cache=True, 93 use_cache=True, 78 94 verbose=True) 79 95 80 print 'create G1' 81 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya') 82 print 'create G2' 83 G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya') 84 print 'create G3' 85 G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') 86 print 'create G4' 87 G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya') 88 print 'add G1+G2+G3+G4' 96 print 'create offshore' 97 G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya') 98 print 'create onshore' 99 G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') 100 print 'create coast' 101 G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya') 102 print 'add' 89 103 G = G1 + G2 + G3 + G4 90 print 'export G'104 print 'export points' 91 105 G.export_points_file(project.combined_dem_name + '.pts') 92 93 #---------------------------------------------------------------------------- ---106 ''' 107 #---------------------------------------------------------------------------- 94 108 # Create the triangular mesh based on overall clipping polygon with a tagged 95 109 # boundary and interior regions defined in project.py along with … … 98 112 99 113 from anuga.pmesh.mesh_interface import create_mesh_from_regions 100 101 region_res = 500000 114 remainder_res = 100000 115 local_res = 15000 116 broome_res = 2500 102 117 coast_res = 500 103 broome_res = 5000 104 interior_regions = [[project.poly_broome, broome_res], 105 [project.poly_region, region_res]] 118 beach_res = 500 119 interior_regions = [[project.poly_broome1, local_res], 120 [project.poly_broome2, broome_res], 121 [project.poly_broome3, coast_res], 122 [project.poly_broome4, beach_res]] 106 123 107 124 print 'number of interior regions', len(interior_regions) 108 109 from anuga.utilities.polygon import plot_polygons 110 if sys.platform == 'win32': 111 #figname = project.outputtimedir + 'pt_hedland_polys' 112 figname = 'broome_polys_test' 113 plot_polygons([project.polyAll,project.poly_broome,project.poly_region], 114 figname, 115 verbose = True) 116 117 print 'start create mesh from regions' 125 print 'estimate of number of triangles', \ 126 number_mesh_triangles(interior_regions, project.polyAll, remainder_res) 127 128 import sys; sys.exit() 129 118 130 from caching import cache 119 131 _ = cache(create_mesh_from_regions, 120 132 project.polyAll, 121 {'boundary_tags': {'topright': [0], 'topleft': [1], 122 'left': [2], 'bottom0': [3], 123 'bottom1': [4], 'bottom2': [5], 124 'bottom3': [6], 'right': [7]}, 125 'maximum_triangle_area': 250000, 126 'filename': meshname, 133 {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2], 134 'e3': [3], 'e4':[4], 'e5': [5], 135 'e6': [6]}, 136 'maximum_triangle_area': remainder_res, 137 'filename': meshname, 127 138 'interior_regions': interior_regions}, 128 verbose = True, evaluate=True) 129 139 verbose = True, evaluate=False) 140 print 'created mesh' 141 import sys; sys.exit() 130 142 #------------------------------------------------------------------------------- 131 143 # Setup computational domain 132 144 #------------------------------------------------------------------------------- 133 domain = Domain(meshname, use_cache = False, verbose = True) 134 135 print domain.statistics() 145 domain = Domain(meshname, use_cache = True, verbose = True) 146 136 147 print 'Number of triangles = ', len(domain) 137 148 print 'The extent is ', domain.get_extent() … … 141 152 domain.set_datadir(project.outputtimedir) 142 153 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) 154 domain.set_minimum_storable_height(0.01) 143 155 144 156 #------------------------------------------------------------------------------- … … 146 158 #------------------------------------------------------------------------------- 147 159 148 tide = 0. 149 #high 150 #tide = 151 #low 152 #tide = 153 160 tide = 0.0 154 161 domain.set_quantity('stage', tide) 155 162 domain.set_quantity('friction', 0.0) 156 print 'hi and file',project.combined_dem_name + '.pts' 157 158 domain.set_quantity('elevation', 159 filename = project.combined_dem_name + '.pts', 163 domain.set_quantity('elevation', 0.0, 164 #filename = project.combined_dem_name + '.pts', 160 165 use_cache = True, 161 166 verbose = True, … … 164 169 165 170 #------------------------------------------------------------------------------- 166 # Setup boundary conditions 167 #------------------------------------------------------------------------------- 168 print 'start ferret2sww'169 # skipped as results in file SU-AU_clipped is correct for all WA 170 171 from anuga. pyvolution.data_manager import ferret2sww171 # Setup boundary conditions 172 #------------------------------------------------------------------------------- 173 ''' 174 print 'start urs2sww' 175 print '', project.boundary_basename 176 from anuga.shallow_water.data_manager import urs2sww 172 177 173 178 south = project.south 174 179 north = project.north 175 west = project.west176 east = project.east180 west = project.west 181 east = project.east 177 182 178 183 #note only need to do when an SWW file for the MOST boundary doesn't exist 179 cache( ferret2sww,184 cache(urs2sww, 180 185 (source_dir + project.boundary_basename, 181 source_dir + project.boundary_basename +'_'+project.basename),186 source_dir + project.boundary_basename), 182 187 {'verbose': True, 183 188 'minlat': south, … … 185 190 'minlon': west, 186 191 'maxlon': east, 187 # 'origin': project.mesh_origin, 188 'origin': domain.geo_reference.get_origin(), 192 #'origin': domain.geo_reference.get_origin(), 189 193 'mean_stage': tide, 190 194 'zscale': 1, #Enhance tsunami 191 195 'fail_on_NaN': False, 192 196 'inverted_bathymetry': True}, 193 197 #evaluate = True, 194 198 verbose = True, 195 dependencies = source_dir + project.boundary_basename + '.sww') 196 199 dependencies = source_dir + project.boundary_basename + '.sww') 200 201 ''' 197 202 print 'Available boundary tags', domain.get_boundary_tags() 198 203 199 Bf = File_boundary(source_dir + project.boundary_basename + '.sww',200 domain, verbose = True)204 #Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 205 # domain, verbose = True) 201 206 Br = Reflective_boundary(domain) 202 207 Bd = Dirichlet_boundary([tide,0,0]) 203 domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left': Bd, 'bottom0': Bd, 204 'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 205 'right': Bd}) 208 209 # 7 min square wave starting at 1 min, 6m high 210 Bw = Time_boundary(domain = domain, 211 f=lambda t: [(60<t<480)*10, 0, 0]) 212 213 domain.set_boundary( {'e0': Bd, 'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd, 214 'e5': Bd, 'e6': Bd} ) 215 206 216 207 217 #------------------------------------------------------------------------------- … … 211 221 t0 = time.time() 212 222 213 for t in domain.evolve(yieldstep = 240, finaltime = 10800):223 for t in domain.evolve(yieldstep = 240, finaltime = 6800): 214 224 domain.write_time() 215 domain.write_boundary_statistics(tags = ' topright')216 217 for t in domain.evolve(yieldstep = 120, finaltime = 16200225 domain.write_boundary_statistics(tags = 'e14') 226 227 for t in domain.evolve(yieldstep = 30, finaltime = 9000 218 228 ,skip_initial_step = True): 219 229 domain.write_time() 220 domain.write_boundary_statistics(tags = ' topright')221 222 for t in domain.evolve(yieldstep = 60, finaltime = 21600230 domain.write_boundary_statistics(tags = 'e14') 231 232 for t in domain.evolve(yieldstep = 240, finaltime = 15000 223 233 ,skip_initial_step = True): 224 234 domain.write_time() 225 domain.write_boundary_statistics(tags = ' topright')235 domain.write_boundary_statistics(tags = 'e14') 226 236 227 for t in domain.evolve(yieldstep = 120, finaltime = 27000228 ,skip_initial_step = True):229 domain.write_time()230 domain.write_boundary_statistics(tags = 'topright')231 232 for t in domain.evolve(yieldstep = 240, finaltime = 36000233 ,skip_initial_step = True):234 domain.write_time()235 domain.write_boundary_statistics(tags = 'topright')236 237 237 print 'That took %.2f seconds' %(time.time()-t0) 238 238 -
anuga_work/production/onslow_2006/project.py
r3769 r3788 5 5 from os import sep, environ, getenv, getcwd 6 6 from os.path import expanduser 7 #from anuga.utilities.polygon import read_polygon 7 from anuga.utilities.polygon import read_polygon, polygon_area 8 8 import sys 9 9 … … 163 163 164 164 polyAll = [d0, d1, d2, d3, d4, d5, d6] 165 165 print 'bounding polygon area', polygon_area(polyAll)/1000000.0 166 166 polygons = [polyAll, export_region] 167 167 figname = 'checking.png' … … 208 208 209 209 poly_onslow = [i0, i1, i2, i3, i4, i5] 210 210 print 'onslow polygon area', polygon_area(poly_onslow)/1000000.0 211 211 #Thevenard Island 212 212 j0 = [294000, 7629000] … … 216 216 217 217 poly_thevenard = [j0, j1, j2, j3] 218 218 print 'thevenard polygon area', polygon_area(poly_thevenard)/1000000.0 219 219 ''' 220 220 # Direction Is … … 237 237 #poly_coast = [l0, l1, l2, l3, l4, l5] 238 238 poly_coast = [l0, l1, l2, l3, l4] 239 239 print 'coast polygon area', polygon_area(poly_coast)/1000000.0 240 240 #general coast and local area to onslow region 241 241 m0 = [270000, 7581000] … … 247 247 248 248 poly_region = [m0, m1, m2, m3, m4, m5] 249 print 'region polygon area', polygon_area(poly_region)/1000000.0 -
anuga_work/production/onslow_2006/run_onslow.py
r3535 r3788 23 23 24 24 # Related major packages 25 from anuga. pyvolution.shallow_water import Domain, Reflective_boundary, \25 from anuga.shallow_water import Domain, Reflective_boundary, \ 26 26 Dirichlet_boundary, Time_boundary, File_boundary 27 from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts 28 from anuga.pyvolution.combine_pts import combine_rectangular_points_files 29 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance 27 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts 28 from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 30 29 from anuga.geospatial_data.geospatial_data import * 31 from anuga. pyvolution.util import Screen_Catcher30 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher 32 31 33 32 # Application specific imports … … 144 143 #------------------------------------------------------------------------------- 145 144 146 domain = pmesh_to_domain_instance(meshname, Domain, 147 use_cache = False, 148 verbose = True) 145 #domain = pmesh_to_domain_instance(meshname, Domain, 146 # use_cache = False, 147 # verbose = True) 148 149 domain = Domain(meshname, use_cache = False, verbose = True) 149 150 150 151 print 'Number of triangles = ', len(domain) -
anuga_work/production/pt_hedland_2006/project.py
r3770 r3788 7 7 #from anuga.utilities.polygon import read_polygon 8 8 import sys 9 from anuga. pmesh.create_mesh import convert_from_latlon_to_utm9 from anuga.coordinate_transforms.redfearn import convert_points_from_latlon_to_utm 10 10 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees 11 11 from time import localtime, strftime -
anuga_work/production/pt_hedland_2006/run_pt_hedland.py
r3535 r3788 19 19 20 20 # Related major packages 21 from anuga. pyvolution.shallow_water import Domain, Reflective_boundary, \21 from anuga.shallow_water import Domain, Reflective_boundary, \ 22 22 Dirichlet_boundary, Time_boundary, File_boundary 23 from anuga. pyvolution.data_manager import convert_dem_from_ascii2netcdf, \23 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \ 24 24 dem2pts 25 from anuga.pyvolution.combine_pts import combine_rectangular_points_files 26 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance 25 from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files 27 26 from shutil import copy 28 27 from os import mkdir, access, F_OK 29 28 from anuga.geospatial_data.geospatial_data import * 30 29 import sys 31 from anuga. pyvolution.util import Screen_Catcher30 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher 32 31 33 32 # Application specific imports -
anuga_work/production/sydney_2006/project.py
r3769 r3788 8 8 import sys 9 9 10 from anuga. pmesh.create_mesh importconvert_from_latlon_to_utm10 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_from_latlon_to_utm 11 11 12 12
Note: See TracChangeset
for help on using the changeset viewer.