Changeset 2615
- Timestamp:
- Mar 28, 2006, 11:00:20 AM (18 years ago)
- Location:
- production/onslow_2006
- Files:
-
- 2 added
- 1 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
production/onslow_2006/project_new.py
r2470 r2615 9 9 10 10 from pmesh.create_mesh import convert_points_from_latlon_to_utm 11 12 from coordinate_transforms.redfearn import degminsec2decimal_degrees 13 14 from time import localtime, strftime 15 16 from os import getcwd 11 17 12 18 #Making assumptions about the location of scenario data 13 19 scenario_dir_name = 'onslow_tsunami_scenario_2006' 14 20 15 16 21 # 250m data to be provided 17 22 coarsename = 'onsl_bathydem250' # get from Neil/Ingo (DEM or topo data) 18 """ 23 19 24 # 30m data to be provided 20 finename = 'onshore_30' # get from Neil/Ingo (DEM or topo data) 21 """ 25 onshore_name = 'onslow_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data) 22 26 23 # clipping region for fine elevation data 24 west = 240000 25 east = 340000 26 south = 7570000 27 north = 7645000 27 offshore_name = 'onslow_offshore_points' 28 29 boundary_basename = 'SU-AU' 28 30 29 31 #swollen/ all data output 30 32 basename = 'source' 33 34 codename = 'project.py' 31 35 32 36 if sys.platform == 'win32': … … 36 40 37 41 #Derive subdirectories and filenames 42 timedir = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir 43 outputdir = home+sep+scenario_dir_name+sep+'output'+sep+timedir+sep 38 44 meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep 39 45 datadir = home+sep+scenario_dir_name+sep+'topographies'+sep 40 outputdir = home+sep+scenario_dir_name+sep+'output'+sep 46 41 47 polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep 42 48 boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep 49 50 codedir = getcwd()+sep 51 52 codedirname = codedir + 'project.py' 43 53 44 54 meshname = meshdir + basename 45 55 46 56 coarsedemname = datadir + coarsename 47 """ 48 finedemname = datadir + finename 49 combineddemname = datadir + 'onslow_combined_elevation' 50 """ 57 58 onshore_dem_name = datadir + onshore_name 59 60 offshore_dem_name = datadir + offshore_name 61 62 combined_dem_name = datadir + 'onslow_combined_elevation' 51 63 52 64 outputname = outputdir + basename #Used by post processing … … 54 66 #!gauge_filename = outputdir + 'onslow_gauges.xya' 55 67 #!gauge_outname = outputdir + 'gauges_max_output.xya' 68 69 # clipping region for fine elevation data 70 eastingmin = 250000 71 eastingmax = 330000 72 northingmin = 7580000 73 northingmax = 7635000 74 75 south = degminsec2decimal_degrees(-22,00,0) 76 north = degminsec2decimal_degrees(-21,10,0) 77 west = degminsec2decimal_degrees(114,30,0) 78 east = degminsec2decimal_degrees(115,30,0) 79 80 # region for visualisation 81 eminviz = 260000 82 emaxviz = 320000 83 nminviz = 7590000 84 nmaxviz = 7630000 56 85 57 86 #Georeferencing … … 97 126 poly_direction = [k0, k1, k2, k3] 98 127 99 #!slump_origin = [385000.0, 6255000.0] #Absolute UTM -
production/onslow_2006/run_onslow_new.py
r2470 r2615 21 21 # Related major packages 22 22 from pyvolution.shallow_water import Domain, Reflective_boundary, \ 23 Dirichlet_boundary, Time_boundary 23 Dirichlet_boundary, Time_boundary, File_boundary 24 24 from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts 25 25 from pyvolution.combine_pts import combine_rectangular_points_files 26 26 from pyvolution.pmesh2domain import pmesh_to_domain_instance 27 from geospatial_data import add_points_files 27 28 28 29 # Application specific imports … … 30 31 from smf import slump_tsunami # Function for submarine mudslide 31 32 33 from shutil import copy 34 from os import mkdir, access, F_OK 35 36 from geospatial_data import * 32 37 33 38 #------------------------------------------------------------------------------- … … 41 46 # filenames 42 47 coarsedemname = project.coarsedemname 43 ''' 44 #finedemname = project.finedemname 45 ''' 48 49 onshore_dem_name = project.onshore_dem_name 50 51 offshore_points = project.offshore_dem_name 52 46 53 meshname = project.meshname+'.msh' 47 54 55 source_dir = project.boundarydir 56 57 # creates copy of code in output dir 58 if access(project.outputdir,F_OK) == 0 : 59 mkdir (project.outputdir) 60 copy (project.codedirname, project.outputdir + project.codename) 61 copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py') 62 63 64 ''' 48 65 # coarse data 49 66 convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True) 50 67 dem2pts(coarsedemname, use_cache=True, verbose=True) 51 ''' 68 69 52 70 # fine data (clipping the points file to smaller area) 53 convert_dem_from_ascii2netcdf( finedemname, use_cache=True, verbose=True)54 dem2pts( finedemname,71 convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True) 72 dem2pts(onshore_dem_name, 55 73 easting_min=project.eastingmin, 56 74 easting_max=project.eastingmax, … … 60 78 verbose=True) 61 79 62 63 # combining the coarse and fine data 64 combine_rectangular_points_files(project.finedemname + '.pts', 65 project.coarsedemname + '.pts', 66 project.combineddemname + '.pts') 67 ''' 68 80 ''' 81 print'create G1' 82 G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya') 83 84 print'create G2' 85 G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts') 86 87 print'add G1+G2' 88 G = G1 + G2 89 90 print'export G' 91 G.new_export_points_file(project.combined_dem_name + '.pts') 92 93 ''' 94 add_points_files( 95 project.offshore_dem_name + '.xya', 96 project.onshore_dem_name + '.pts', 97 project.combined_dem_name + '.pts') 98 ''' 69 99 #------------------------------------------------------------------------------- 70 100 # Create the triangular mesh based on overall clipping polygon with a tagged … … 73 103 #------------------------------------------------------------------------------- 74 104 75 from pmesh. create_meshimport create_mesh_from_regions105 from pmesh.mesh_interface import create_mesh_from_regions 76 106 77 107 # original … … 80 110 [project.poly_thevenard, interior_res], 81 111 [project.poly_direction, interior_res]] 82 83 #FIXME: Fix caching of this one once the mesh_interface is ready 112 #[project.testpoly, interior_res]] 113 print 'number of interior regions', len(interior_regions) 114 84 115 from caching import cache 85 116 _ = cache(create_mesh_from_regions, … … 88 119 'left': [2], 'bottom': [3], 89 120 'bottomright': [4], 'topright': [5]}, 90 ' resolution': 100000,121 'maximum_triangle_area': 100000, 91 122 'filename': meshname, 92 'interior_regions': interior_regions, 93 'UTM': True, 94 'refzone': project.refzone}, 123 'interior_regions': interior_regions}, 95 124 verbose = True) 96 125 … … 106 135 print 'Number of triangles = ', len(domain) 107 136 print 'The extent is ', domain.get_extent() 137 print domain.statistics() 108 138 109 139 domain.set_name(project.basename) 110 140 domain.set_datadir(project.outputdir) 111 domain.set_quantities_to_be_stored(['stage' , 'xmomentum', 'ymomentum'])141 domain.set_quantities_to_be_stored(['stage']) 112 142 113 143 … … 136 166 domain.set_quantity('stage', tide) 137 167 domain.set_quantity('friction', 0.0) 168 print 'hi and file',project.combined_dem_name + '.pts' 138 169 domain.set_quantity('elevation', 139 # filename = project.combineddemname + '.pts', 140 filename = project.coarsedemname + '.pts', 141 use_cache = True, 142 verbose = True 170 # 0. 171 # filename = project.onshore_dem_name + '.pts', 172 filename = project.combined_dem_name + '.pts', 173 # filename = project.offshore_dem_name + '.pts', 174 use_cache = False, 175 verbose = True, 176 alpha = 0.1 143 177 ) 144 178 print 'hi1' 145 179 146 180 #------------------------------------------------------------------------------- … … 156 190 157 191 cache(ferret2sww, 158 (source_dir +project.boundary_basename,159 project.boundary_basename),192 (source_dir + project.boundary_basename, 193 source_dir + project.boundary_basename), 160 194 {'verbose': True, 161 'minlat': south-1, 162 'maxlat': north+1, 163 'minlon': west-1, 164 'maxlon': east+1, 165 'origin': project.mesh_origin, 195 # note didn't work with the below 196 # 'minlat': south - 1, 197 # 'maxlat': north + 1, 198 # 'minlon': west - 1, 199 # 'maxlon': east + 1, 200 'minlat': south, 201 'maxlat': north, 202 'minlon': west, 203 'maxlon': east, 204 # 'origin': project.mesh_origin, 205 'origin': domain.geo_reference.get_origin(), 166 206 'mean_stage': tide, 167 207 'zscale': 1, #Enhance tsunami … … 174 214 print 'Available boundary tags', domain.get_boundary_tags() 175 215 216 Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 217 domain, verbose = True) 176 218 Br = Reflective_boundary(domain) 177 219 Bd = Dirichlet_boundary([tide,0,0]) … … 182 224 f=lambda t: [(60<t<480)*6, 0, 0]) 183 225 184 domain.set_boundary( {'top': B w, 'topleft': Br,226 domain.set_boundary( {'top': Bf, 'topleft': Bf, 185 227 'left': Br, 'bottom': Br, 186 'bottomright': Br, 'topright': B r} )228 'bottomright': Br, 'topright': Bf} ) 187 229 188 230 … … 190 232 # Evolve system through time 191 233 #------------------------------------------------------------------------------- 192 193 234 import time 194 235 t0 = time.time() … … 199 240 200 241 print 'That took %.2f seconds' %(time.time()-t0) 242 243 print 'finished'
Note: See TracChangeset
for help on using the changeset viewer.