Changeset 9321
- Timestamp:
- Sep 4, 2014, 11:41:04 AM (11 years ago)
- Location:
- trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts
- Files:
-
- 3 edited
-
project.py (modified) (4 diffs)
-
raster_outputs.py (modified) (3 diffs)
-
run_model.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py
r9320 r9321 34 34 ## Total model runtime in seconds (starts at time=0.) 35 35 # 36 finalTime= 4.0*3600.036 finalTime=0.01*3600.0 37 37 38 38 ## Model time between output file timesteps in seconds … … 137 137 ######################################################################### 138 138 139 ## Output tif details 140 # 141 # pixle size (m) for output tifs 142 output_tif_cellsize = 250. 143 # EPSG CODE (projection information) for outputs tifs 144 EPSG_CODE = 32647 145 139 146 ## Storage of max quantities 140 147 # … … 153 160 # before storage, leading to smaller file sizes but a less accurate 154 161 # representation of what the algorithm is doing 155 store_vertices_uniquely= True162 store_vertices_uniquely=False 156 163 157 164 ## Checkpointing … … 159 166 # Need to restructure to allow checkpointing 160 167 # checkpoint_time=30 161 162 168 163 169 #@@@@@@@@@@@@@@@@@# -
trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/raster_outputs.py
r9320 r9321 1 1 """ 2 Code to make spatial outputs 2 3 Code to make spatial outputs (GeoTifs) 3 4 4 5 """ … … 11 12 from anuga.utilities import plot_utils as util 12 13 from anuga.utilities import spatialInputUtil as su 13 import project # Get the bounding polygon14 15 16 #@@@@@@@@@@@@@@@@@@@@#17 # START EDITING HERE #18 #@@@@@@@@@@@@@@@@@@@@#19 20 ###########################################################21 ## INPUT DATA22 23 #@ Name of sww to read outputs from24 swwFile='./MODEL_OUTPUTS/RUN_20140902_205219_Ondoy_firstOrder_150m/Ondoy_firstOrder_150m.sww'25 26 #@ myTimeStep option from util.Make_Geotif27 # use 'max' to plot the maxima28 # use [0, 5, 10] to plot the29 # use 'collected_max' to read the csv outputs from collect_max_quantities_operator30 myTimeStep='collected_max'31 #myTimeStep='max'32 33 #@ Write outputs to this folder (INCLUDE TRAILING SLASH /)34 tif_output_dir=os.path.dirname(swwFile)+'/TIFS/'35 36 #@ Desired raster cellSize (m)37 CellSize=5.038 39 #@ EPSG projection info40 EPSG_CODE=312341 42 #@ Elevation raster for 'high-res-drape' depth plot43 elevation_raster='./../../../../../DATA/Manila_DEM_20121008/DEM_Mosaic/Processing/new_tiffs/new_mosaic.vrt'44 #elevation_raster='../../../../../DATA/Manila_DSM/manila_DSM/Processed/manila_DSM.vrt'45 46 #@ Depth threshold for high-res-drape' plot (m). If ANUGA's triangle has depth47 # < depth_threshold, then depth=0 in all high-res cells in that triangle48 depth_threshold=0.0549 50 #@ Polygon to clip 'high-res-drape' plot. Must be provided (e.g. can use bounding polygon or another choice)51 clipPolygon='../../GIS/PlotClipPolygon/PlotClipPolygon.shp'52 53 #@ Layer for above (usually shapefile name without .shp)54 clipPolygonLayer='PlotClipPolygon'55 56 #@ Skip the creation of rasters from sww or collected-max-quantity files.57 # Useful for debugging, or changing the DEM (but by default old files are overwritten).58 skipInitialRasterCreation=False59 14 60 15 #@@@@@@@@@@@@@@@@@@# 61 16 # STOP EDITING HERE 62 17 #@@@@@@@@@@@@@@@@@@# 63 64 ############################################################ 65 66 # Make the geotifs 67 if(not skipInitialRasterCreation): 68 if(myTimeStep=='collected_max'): 69 # Read max quantity output files inputs 70 maxQFiles=glob.glob(os.path.dirname(swwFile)+'/*_UH_MAX.csv') 71 if(len(maxQFiles)==0): 72 raise Exception, 'Cannot find any files containing collected maxima' 73 for i in range(len(maxQFiles)): 74 if(i==0): 75 maxQs=numpy.genfromtxt(maxQFiles[i], delimiter=',') 76 else: 77 maxQs=numpy.vstack([maxQs, numpy.genfromtxt(maxQFiles[i], delimiter=',')]) 78 # Make the geotiff's 79 for i, quant in enumerate(['stage', 'depth', 'velocity', 'depthIntegratedVelocity']): 80 # FIXME: The interpolation function is remade for every quantity, since 81 # only a 3 column array can be passed to Make_Geotif 82 # Could make it fast (do it only once) by changing Make_Geotif 83 tmpArr=maxQs[:,[0, 1, i+2]] 84 util.Make_Geotif(tmpArr, output_quantities=[quant+'_MAX'], 85 CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 86 bounding_polygon=project.bounding_polygon, output_dir=tif_output_dir) 87 else: 88 util.Make_Geotif(swwFile, myTimeStep=myTimeStep, 89 output_quantities=['depth', 'stage', 'elevation', 'velocity', 'depthIntegratedVelocity', 'friction'], 90 CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, bounding_polygon=project.bounding_polygon, 91 output_dir=tif_output_dir) 18 # 19 # 92 20 93 21 def make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir): … … 157 85 return 158 86 87 ############################################################ 159 88 160 # Get extent of geotifs 161 sampleRast=glob.glob(tif_output_dir+'*.tif')[0] 162 rasterExtent=su.getRasterExtent(sampleRast) 89 def make_me_some_tifs(swwFile, 90 bounding_polygon, 91 EPSG_CODE, 92 myTimeStep='collected_max', 93 tif_output_subdir='/TIFS/', 94 CellSize=5.0, 95 make_highres_drape_plot=False, 96 elevation_raster=None, 97 depth_threshold=None, 98 clipPolygon=None, 99 clipPolygonLayer=None 100 ): 101 """ 163 102 164 # Make the resampled elevation data 165 make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir) 103 ### INPUT DATA 104 - swwFile -- Full path name of sww to read outputs from 105 - bounding_polygon -- ANUGA's bounding polygon, or another polygon to clip the rasters to 106 - EPSG_CODE -- EPSG projection info as an integer like 3123 [which is over Manila, Philippines] 107 - myTimeStep -- option from util.Make_Geotif 108 use 'max' to plot the maxima 109 use [0, 5, 10] to plot the 110 use 'collected_max' to read the csv outputs from collect_max_quantities_operator 111 - tif_outputdir -- Write outputs to this folder inside the swwFile's directory (MUST INCLUDE TRAILING SLASH /) 112 - CellSize -- Desired raster cellSize (m) 113 - make_highres_drape_plot -- True/False, Make a high-res drape plot? 114 - elevation_raster -- Filename of elevation raster for 'high-res-drape' 115 depth plot [from subtracting stage from high res topography] 116 - depth_threshold -- Depth threshold for high-res-drape' plot (m). If ANUGA's triangle has depth 117 < depth_threshold, then depth=0 in all high-res cells in that triangle 118 - clipPolygon -- Polygon to clip 'high-res-drape' plot. Must be 119 provided if make_highres_drape_plot==True (can use bounding polygon or another choice) 120 - clipPolygonLayer -- Layer for above as used by gdal (usually shapefile name without .shp) 121 122 ## OUTPUT 123 Nothing is returned, but tifs are made in tif_output_subdir inside the swwFile directory 166 124 167 elevation=glob.glob(tif_output_dir+'LIDAR_resampled*')[0] 168 mask=glob.glob(tif_output_dir+'Mask_resampled*')[0] 169 if(myTimeStep=='collected_max'): 170 depth=glob.glob(tif_output_dir+'PointData_depth_*')[0] 171 stage=glob.glob(tif_output_dir+'PointData_stage_*')[0] 172 else: 173 depth=glob.glob(tif_output_dir+'*_depth_max.tif')[0] 174 stage=glob.glob(tif_output_dir+'*_stage_max.tif')[0] 125 """ 175 126 176 # Call gdal_calc 177 gdal_calc_command(stage, depth, elevation, mask, depth_threshold) 127 tif_output_dir=os.path.dirname(swwFile)+tif_output_subdir 128 # Make the geotifs 129 if(myTimeStep=='collected_max'): 130 # Read max quantity output files inputs 131 maxQFiles=glob.glob(os.path.dirname(swwFile)+'/*_UH_MAX.csv') 132 if(len(maxQFiles)==0): 133 raise Exception, 'Cannot find any files containing collected maxima' 134 for i in range(len(maxQFiles)): 135 if(i==0): 136 maxQs=numpy.genfromtxt(maxQFiles[i], delimiter=',') 137 else: 138 maxQs=numpy.vstack([maxQs, numpy.genfromtxt(maxQFiles[i], delimiter=',')]) 139 # Make the geotiff's 140 for i, quant in enumerate(['stage', 'depth', 'velocity', 'depthIntegratedVelocity']): 141 # FIXME: The interpolation function is remade for every quantity, since 142 # only a 3 column array can be passed to Make_Geotif 143 # Could make it fast (do it only once) by changing Make_Geotif 144 tmpArr=maxQs[:,[0, 1, i+2]] 145 util.Make_Geotif(tmpArr, output_quantities=[quant+'_MAX'], 146 CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 147 bounding_polygon=bounding_polygon, output_dir=tif_output_dir) 148 149 # Also plot elevation + friction 150 util.Make_Geotif(swwFile, myTimeStep=[0], 151 output_quantities=['elevation', 'friction'], 152 CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 153 bounding_polygon=bounding_polygon, 154 output_dir=tif_output_dir) 155 else: 156 util.Make_Geotif(swwFile, myTimeStep=myTimeStep, 157 output_quantities=['depth', 'stage', 'elevation', 'velocity', 'depthIntegratedVelocity', 'friction'], 158 CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 159 bounding_polygon=bounding_polygon, 160 output_dir=tif_output_dir) 161 162 # Early finish 163 if not make_highres_drape_plot: 164 return 165 166 # Get extent of geotifs 167 sampleRast=glob.glob(tif_output_dir+'*.tif')[0] 168 rasterExtent=su.getRasterExtent(sampleRast) 169 170 # Make the resampled elevation data 171 make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir) 172 173 elevation=glob.glob(tif_output_dir+'LIDAR_resampled*')[0] 174 mask=glob.glob(tif_output_dir+'Mask_resampled*')[0] 175 if(myTimeStep=='collected_max'): 176 depth=glob.glob(tif_output_dir+'PointData_depth_*')[0] 177 stage=glob.glob(tif_output_dir+'PointData_stage_*')[0] 178 else: 179 depth=glob.glob(tif_output_dir+'*_depth_max.tif')[0] 180 stage=glob.glob(tif_output_dir+'*_stage_max.tif')[0] 181 182 # Call gdal_calc 183 gdal_calc_command(stage, depth, elevation, mask, depth_threshold) 184 185 return -
trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py
r9320 r9321 31 31 # Routines defined by the user 32 32 import user_functions 33 34 # Functions to make geotifs 35 import raster_outputs 33 36 34 37 # Record the time so we can report how long the simulation takes … … 120 123 121 124 # Merge the parallel SWW files together 125 os.chdir(project.output_dir) 122 126 if myid == 0 and numprocs>1: 123 127 print 'Number of processors %g ' %numprocs … … 127 131 print 'Broadcast time %.2f seconds'%domain.communication_broadcast_time 128 132 129 os.chdir(project.output_dir)130 133 anuga.utilities.sww_merge.sww_merge_parallel(project.scenario, np=numprocs, verbose=True, delete_old=True) 131 134 132 135 # FIXME: Make geotif rasters 136 if myid == 0: 137 try: 138 raster_outputs.make_me_some_tifs(swwFile = './'+project.scenario+'.sww', 139 bounding_polygon=project.bounding_polygon, 140 EPSG_CODE = project.EPSG_CODE, 141 CellSize= project.output_tif_cellsize) 142 except: 143 print 'GeoTif creation failed, you can try manually using raster_outputs.py or anuga.utilities.plot_utils.Make_Geotif' 144 133 145 146 barrier() 134 147 finalize() 135 148
Note: See TracChangeset
for help on using the changeset viewer.
