[5005] | 1 | """Script for running tsunami inundation scenario for Dampier, WA, Australia. |
---|
| 2 | |
---|
| 3 | Source data such as elevation and boundary data is assumed to be available in |
---|
| 4 | directories specified by project.py |
---|
| 5 | The output sww file is stored in project.output_time_dir |
---|
| 6 | |
---|
| 7 | The scenario is defined by a triangular mesh created from project.polygon, |
---|
| 8 | the elevation data and a simulated submarine landslide. |
---|
| 9 | |
---|
| 10 | Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006 |
---|
| 11 | """ |
---|
| 12 | |
---|
| 13 | #------------------------------------------------------------------------------ |
---|
| 14 | # Import necessary modules |
---|
| 15 | #------------------------------------------------------------------------------ |
---|
| 16 | |
---|
| 17 | # Standard modules |
---|
| 18 | from os import sep |
---|
| 19 | from os.path import dirname, basename |
---|
| 20 | from os import mkdir, access, F_OK |
---|
| 21 | from shutil import copy |
---|
| 22 | import time |
---|
| 23 | import sys |
---|
| 24 | |
---|
| 25 | # Related major packages |
---|
| 26 | from anuga.shallow_water import Domain |
---|
[5012] | 27 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,start_screen_catcher, copy_code_files |
---|
[5005] | 28 | from anuga.geospatial_data.geospatial_data import * |
---|
| 29 | |
---|
| 30 | # Application specific imports |
---|
| 31 | import project # Definition of file names and polygons |
---|
| 32 | |
---|
| 33 | #------------------------------------------------------------------------------ |
---|
| 34 | # Copy scripts to time stamped output directory and capture screen |
---|
| 35 | # output to file |
---|
| 36 | #------------------------------------------------------------------------------ |
---|
| 37 | |
---|
[5012] | 38 | print 'time stamp: ',project.build_time |
---|
[5005] | 39 | |
---|
| 40 | copy_code_files(project.output_build_time_dir,__file__, |
---|
| 41 | dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
| 42 | |
---|
[5012] | 43 | start_screen_catcher(project.output_build_time_dir) |
---|
[5005] | 44 | |
---|
[5150] | 45 | print 'time stamp: ',project.build_time |
---|
| 46 | print 'USER: ', project.user, project.host |
---|
[5005] | 47 | |
---|
[5150] | 48 | |
---|
[5005] | 49 | #------------------------------------------------------------------------------- |
---|
| 50 | # Preparation of topographic data |
---|
| 51 | # |
---|
| 52 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 53 | # Do for coarse and fine data |
---|
| 54 | # Fine pts file to be clipped to area of interest |
---|
| 55 | #------------------------------------------------------------------------------- |
---|
| 56 | |
---|
| 57 | # topography directory filenames |
---|
| 58 | onshore_in_dir_name = project.onshore_in_dir_name |
---|
| 59 | coast_in_dir_name = project.coast_in_dir_name |
---|
[5012] | 60 | island_in_dir_name = project.island_in_dir_name |
---|
[5005] | 61 | offshore_in_dir_name = project.offshore_in_dir_name |
---|
| 62 | |
---|
| 63 | onshore_dir_name = project.onshore_dir_name |
---|
| 64 | coast_dir_name = project.coast_dir_name |
---|
[5012] | 65 | island_dir_name = project.island_dir_name |
---|
[5005] | 66 | offshore_dir_name = project.offshore_dir_name |
---|
| 67 | |
---|
| 68 | # creates DEM from asc data |
---|
| 69 | print "creates DEMs from ascii data" |
---|
| 70 | convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True) |
---|
[5012] | 71 | convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True) |
---|
[5005] | 72 | |
---|
| 73 | #creates pts file for onshore DEM |
---|
| 74 | print "creates pts file for onshore DEM" |
---|
| 75 | dem2pts(onshore_dir_name, |
---|
| 76 | use_cache=True, |
---|
| 77 | verbose=True) |
---|
| 78 | |
---|
| 79 | #creates pts file for island DEM |
---|
[5012] | 80 | dem2pts(island_dir_name, use_cache=True, verbose=True) |
---|
[5005] | 81 | |
---|
| 82 | print'create Geospatial data1 objects from topographies' |
---|
[5150] | 83 | G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True) |
---|
[5012] | 84 | print'finished reading in %s.pts' %onshore_dir_name |
---|
[5150] | 85 | G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True) |
---|
| 86 | G3 = Geospatial_data(file_name = island_dir_name + '.pts',verbose=True) |
---|
| 87 | G4 = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True) |
---|
[5012] | 88 | print'finished reading in files' |
---|
[5005] | 89 | |
---|
[5150] | 90 | #G_onshore_clip = G1.clip(project.poly_all, verbose=True) |
---|
| 91 | #print 'finished clip' |
---|
| 92 | #G_onshore_clip_clipout=G_onshore_clip.clip_outside(project.poly_cbd, verbose=True) |
---|
| 93 | #print 'finished clipout' |
---|
| 94 | ## reduce resolution by 10 time (1/10) of the original file. This was determined |
---|
| 95 | ## as the resolution for this region is going to be 20000 and the asc grid cellsize is |
---|
| 96 | ## 10m there (100). There is 2 20000m triangles in a 200m x 200m grid sqrt of 40000 is about 200, |
---|
| 97 | ## a 200m x 200m area would contain 400 points at a 10m grid spacing and we only need 4. So 4/400 is 1/100 or 0.01 |
---|
| 98 | ## 1/100 reduction is possible but a 1/10 is enough to help with file size and loading |
---|
| 99 | #G_onshore_clip_clipout_reduced=G_onshore_clip_clipout.split(0.1, verbose=True) |
---|
| 100 | # |
---|
| 101 | #G_cbd = G1.clip(projec.poly_cbd,verbose=True) |
---|
| 102 | ##G_cbd_reduced = G_cbd.split(0.25, verbose=True) |
---|
| 103 | # |
---|
| 104 | #G_bathy = G4.clip(project.poly_all) |
---|
| 105 | #G_island = G3.clip(project.poly_all) |
---|
| 106 | #G_coast = G2.clip(project.poly_all) |
---|
| 107 | # |
---|
| 108 | #print'add all geospatial objects' |
---|
| 109 | #G = G_onshore_clip_clipout_reduced + G_cbd_reduced + G_bathy + G_island + G_coast |
---|
[5012] | 110 | |
---|
[5150] | 111 | G = G1 + G2 + G3 + G4 |
---|
[5012] | 112 | |
---|
[5150] | 113 | G_clip = G.clip(project.poly_all, verbose=True) |
---|
[5012] | 114 | |
---|
[5005] | 115 | print'clip combined geospatial object by bounding polygon' |
---|
| 116 | #G_clipped = G.clip(project.bounding_polygon) |
---|
| 117 | #FIXME: add a clip function to pts |
---|
| 118 | #print'shape of clipped data', G_clipped.get_data_points().shape |
---|
| 119 | |
---|
| 120 | print'export combined DEM file' |
---|
| 121 | if access(project.topographies_dir,F_OK) == 0: |
---|
| 122 | mkdir (project.topographies_dir) |
---|
[5150] | 123 | G_clip.export_points_file(project.combined_dir_name + '.txt') |
---|
[5005] | 124 | #G_clipped.export_points_file(project.combined_dir_name + '.xya') |
---|
[5150] | 125 | G_small,G_other = G_clip.split(0.01, verbose=True) |
---|
[5012] | 126 | G_small.export_points_file(project.combined_dir_name_small + '.txt') |
---|
[5005] | 127 | |
---|
| 128 | ''' |
---|
| 129 | print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya' |
---|
| 130 | G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya') |
---|
| 131 | print'split' |
---|
| 132 | G_all_1, G_all_2 = G_all.split(.10) |
---|
| 133 | print'export 1' |
---|
| 134 | G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya') |
---|
| 135 | print'export 2' |
---|
| 136 | G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya') |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | #------------------------------------------------------------------------- |
---|
| 140 | # Convert URS to SWW file for boundary conditions |
---|
| 141 | #------------------------------------------------------------------------- |
---|
| 142 | print 'starting to create boundary conditions' |
---|
| 143 | boundaries_in_dir_name = project.boundaries_in_dir_name |
---|
| 144 | |
---|
| 145 | from anuga.shallow_water.data_manager import urs2sww |
---|
| 146 | |
---|
| 147 | print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary |
---|
| 148 | print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary |
---|
| 149 | |
---|
| 150 | #import sys; sys.exit() |
---|
| 151 | |
---|
| 152 | #if access(project.boundaries_dir,F_OK) == 0: |
---|
| 153 | # mkdir (project.boundaries_dir) |
---|
| 154 | |
---|
| 155 | from caching import cache |
---|
| 156 | cache(urs2sww, |
---|
| 157 | (boundaries_in_dir_name, |
---|
| 158 | project.boundaries_dir_name1), |
---|
| 159 | {'verbose': True, |
---|
| 160 | 'minlat': project.south_boundary, |
---|
| 161 | 'maxlat': project.north_boundary, |
---|
| 162 | 'minlon': project.west_boundary, |
---|
| 163 | 'maxlon': project.east_boundary, |
---|
| 164 | # 'minlat': project.south, |
---|
| 165 | # 'maxlat': project.north, |
---|
| 166 | # 'minlon': project.west, |
---|
| 167 | # 'maxlon': project.east, |
---|
| 168 | 'mint': 0, 'maxt': 40000, |
---|
| 169 | # 'origin': domain.geo_reference.get_origin(), |
---|
| 170 | 'mean_stage': project.tide, |
---|
| 171 | # 'zscale': 1, #Enhance tsunami |
---|
| 172 | 'fail_on_NaN': False}, |
---|
| 173 | verbose = True, |
---|
| 174 | ) |
---|
| 175 | # dependencies = source_dir + project.boundary_basename + '.sww') |
---|
| 176 | |
---|
| 177 | ''' |
---|
| 178 | |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | |
---|
| 184 | |
---|