[3827] | 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 | |
---|
| 26 | # Related major packages |
---|
| 27 | from anuga.shallow_water import Domain |
---|
| 28 | from anuga.shallow_water import Dirichlet_boundary |
---|
| 29 | from anuga.shallow_water import File_boundary |
---|
| 30 | from anuga.shallow_water import Reflective_boundary |
---|
| 31 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
| 32 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 33 | from anuga.geospatial_data.geospatial_data import * |
---|
[3870] | 34 | from anuga.abstract_2d_finite_volumes.util import Screen_Catcher |
---|
[3827] | 35 | |
---|
| 36 | # Application specific imports |
---|
| 37 | import project # Definition of file names and polygons |
---|
| 38 | |
---|
| 39 | #------------------------------------------------------------------------------ |
---|
| 40 | # Copy scripts to time stamped output directory and capture screen |
---|
| 41 | # output to file |
---|
| 42 | #------------------------------------------------------------------------------ |
---|
[3853] | 43 | |
---|
[3851] | 44 | if access(project.output_build_time_dir,F_OK) == 0: |
---|
| 45 | mkdir (project.output_build_time_dir) |
---|
[3827] | 46 | copy (dirname(project.__file__) +sep+ project.__name__+'.py', |
---|
[3851] | 47 | project.output_build_time_dir + project.__name__+'.py') #copies project.py |
---|
| 48 | copy (__file__, project.output_build_time_dir + basename(__file__)) |
---|
| 49 | print 'files '+ project.__name__+'.py and '+ basename(__file__)+' copied to '+ project.output_build_time_dir #copies this file |
---|
| 50 | #import sys; sys.exit() |
---|
[3827] | 51 | |
---|
[3853] | 52 | |
---|
| 53 | #normal screen output is stored in |
---|
[3871] | 54 | screen_output_name = project.output_build_time_dir + "screen_output.txt" |
---|
| 55 | screen_error_name = project.output_build_time_dir + "screen_error.txt" |
---|
[3853] | 56 | |
---|
| 57 | #used to catch screen output to file |
---|
| 58 | sys.stdout = Screen_Catcher(screen_output_name) |
---|
| 59 | sys.stderr = Screen_Catcher(screen_error_name) |
---|
| 60 | |
---|
| 61 | print 'USER: ', project.user |
---|
| 62 | |
---|
[3827] | 63 | #------------------------------------------------------------------------------- |
---|
| 64 | # Preparation of topographic data |
---|
| 65 | # |
---|
| 66 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 67 | # Do for coarse and fine data |
---|
| 68 | # Fine pts file to be clipped to area of interest |
---|
| 69 | #------------------------------------------------------------------------------- |
---|
| 70 | |
---|
| 71 | # topography directory filenames |
---|
| 72 | onshore_dir_name = project.onshore_dir_name |
---|
| 73 | coast_dir_name = project.coast_dir_name |
---|
| 74 | islands_dir_name = project.islands_dir_name |
---|
| 75 | offshore_dir_name = project.offshore_dir_name |
---|
| 76 | offshore_dir_name1 = project.offshore_dir_name1 |
---|
| 77 | offshore_dir_name2 = project.offshore_dir_name2 |
---|
| 78 | offshore_dir_name3 = project.offshore_dir_name3 |
---|
| 79 | offshore_dir_name4 = project.offshore_dir_name4 |
---|
| 80 | offshore_dir_name5 = project.offshore_dir_name5 |
---|
| 81 | offshore_dir_name6 = project.offshore_dir_name6 |
---|
| 82 | offshore_dir_name7 = project.offshore_dir_name7 |
---|
| 83 | offshore_dir_name8 = project.offshore_dir_name8 |
---|
| 84 | offshore_dir_name9 = project.offshore_dir_name9 |
---|
| 85 | offshore_dir_name10 = project.offshore_dir_name10 |
---|
| 86 | offshore_dir_name11 = project.offshore_dir_name11 |
---|
| 87 | offshore_dir_name12 = project.offshore_dir_name12 |
---|
| 88 | offshore_dir_name13 = project.offshore_dir_name13 |
---|
| 89 | offshore_dir_name14 = project.offshore_dir_name14 |
---|
| 90 | |
---|
| 91 | # creates DEM from asc data |
---|
| 92 | convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True) |
---|
| 93 | convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True) |
---|
| 94 | |
---|
| 95 | #creates pts file for onshore DEM |
---|
| 96 | dem2pts(onshore_dir_name, |
---|
| 97 | # easting_min=project.eastingmin, |
---|
| 98 | # easting_max=project.eastingmax, |
---|
| 99 | # northing_min=project.northingmin, |
---|
| 100 | # northing_max= project.northingmax, |
---|
| 101 | use_cache=True, |
---|
| 102 | verbose=True) |
---|
| 103 | |
---|
| 104 | #creates pts file for islands DEM |
---|
| 105 | dem2pts(islands_dir_name, use_cache=True, verbose=True) |
---|
| 106 | |
---|
[3851] | 107 | print'create Geospatial data objects from topographies' |
---|
[3827] | 108 | G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts') |
---|
| 109 | G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya') |
---|
| 110 | G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts') |
---|
| 111 | G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya') |
---|
| 112 | G_off1 = Geospatial_data(file_name = project.offshore_dir_name1 + '.xya') |
---|
| 113 | G_off2 = Geospatial_data(file_name = project.offshore_dir_name2 + '.xya') |
---|
| 114 | G_off3 = Geospatial_data(file_name = project.offshore_dir_name3 + '.xya') |
---|
| 115 | G_off4 = Geospatial_data(file_name = project.offshore_dir_name4 + '.xya') |
---|
| 116 | G_off5 = Geospatial_data(file_name = project.offshore_dir_name5 + '.xya') |
---|
| 117 | G_off6 = Geospatial_data(file_name = project.offshore_dir_name6 + '.xya') |
---|
| 118 | G_off7 = Geospatial_data(file_name = project.offshore_dir_name7 + '.xya') |
---|
| 119 | G_off8 = Geospatial_data(file_name = project.offshore_dir_name8 + '.xya') |
---|
| 120 | G_off9 = Geospatial_data(file_name = project.offshore_dir_name9 + '.xya') |
---|
| 121 | G_off10 = Geospatial_data(file_name = project.offshore_dir_name10 + '.xya') |
---|
| 122 | G_off11 = Geospatial_data(file_name = project.offshore_dir_name11 + '.xya') |
---|
| 123 | G_off12 = Geospatial_data(file_name = project.offshore_dir_name12 + '.xya') |
---|
| 124 | G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya') |
---|
| 125 | G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya') |
---|
| 126 | |
---|
[3851] | 127 | print'add all geospatial objects' |
---|
[3827] | 128 | G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \ |
---|
| 129 | + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \ |
---|
| 130 | + G_off13 + G_off14 |
---|
| 131 | |
---|
[3851] | 132 | print'clip combined geospatial object by bounding polygon' |
---|
[3827] | 133 | G.clip(project.bounding_polygon) |
---|
| 134 | #FIXME: add a clip function to pts |
---|
| 135 | |
---|
[3851] | 136 | print'export combined DEM file' |
---|
| 137 | if access(project.topographies_time_dir,F_OK) == 0: |
---|
| 138 | mkdir (project.topographies_time_dir) |
---|
| 139 | G.export_points_file(project.combined_time_dir_name + '.pts') |
---|
[3827] | 140 | |
---|
| 141 | #------------------------------------------------------------------------- |
---|
| 142 | # Convert URS to SWW file for boundary conditions |
---|
| 143 | #------------------------------------------------------------------------- |
---|
[3851] | 144 | print 'starting to create boundary conditions' |
---|
[3853] | 145 | boundaries_in_dir_name = project.boundaries_in_dir_name |
---|
[3827] | 146 | |
---|
[3851] | 147 | from anuga.shallow_water.data_manager import urs2sww |
---|
[3827] | 148 | |
---|
[3851] | 149 | print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary |
---|
| 150 | print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary |
---|
[3827] | 151 | |
---|
[3851] | 152 | if access(project.boundaries_time_dir,F_OK) == 0: |
---|
| 153 | mkdir (project.boundaries_time_dir) |
---|
[3877] | 154 | urs2sww(boundaries_in_dir_name,basename_out= project.boundaries_time_dir_name, |
---|
| 155 | minlat=project.south_boundary, maxlat=project.north_boundary, |
---|
| 156 | minlon= project.west_boundary, maxlon=project.east_boundary, |
---|
| 157 | mint=0, maxt= 35000, |
---|
| 158 | verbose='true') |
---|
| 159 | ''' |
---|
| 160 | from caching import cache |
---|
[3871] | 161 | cache(urs2sww, |
---|
| 162 | (boundaries_in_dir_name, |
---|
[3877] | 163 | project.boundaries_time_dir_name), |
---|
[3871] | 164 | {'verbose': True, |
---|
| 165 | 'minlat': project.south_boundary, |
---|
| 166 | 'maxlat': project.north_boundary, |
---|
| 167 | 'minlon': project.west_boundary, |
---|
| 168 | 'maxlon': project.east_boundary, |
---|
| 169 | 'mint': 0, 'maxt': 35000, |
---|
[3877] | 170 | # 'origin': domain.geo_reference.get_origin(), |
---|
[3871] | 171 | 'mean_stage': project.tide, |
---|
| 172 | # 'zscale': 1, #Enhance tsunami |
---|
[3877] | 173 | 'fail_on_NaN': False}, |
---|
[3871] | 174 | verbose = True, |
---|
| 175 | ) |
---|
| 176 | # dependencies = source_dir + project.boundary_basename + '.sww') |
---|
[3877] | 177 | ''' |
---|
[3827] | 178 | |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | |
---|
[3851] | 184 | |
---|
| 185 | |
---|