[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 * |
---|
| 34 | |
---|
| 35 | # Application specific imports |
---|
| 36 | import project # Definition of file names and polygons |
---|
| 37 | |
---|
| 38 | #------------------------------------------------------------------------------ |
---|
| 39 | # Copy scripts to time stamped output directory and capture screen |
---|
| 40 | # output to file |
---|
| 41 | #------------------------------------------------------------------------------ |
---|
[3851] | 42 | #print 'project.build_time', project.build_time |
---|
| 43 | #print 'project.output_build_time_dir',project.output_build_time_dir |
---|
| 44 | #print '__file__, project.output_time_dir + basename(__file__)', __file__, project.output_build_time_dir + basename(__file__) |
---|
| 45 | #print 'dirname(project.__file__) +sep+ project.__name__', dirname(project.__file__) +sep+ project.__name__+'.py' |
---|
| 46 | #print'project.topographies_time_dir',project.topographies_time_dir, project.combined_time_dir_name + '.pts' |
---|
| 47 | #print project.boundaries_time_dir, project.boundaries_time_dir_name |
---|
| 48 | #import sys; sys.exit() |
---|
| 49 | if access(project.output_build_time_dir,F_OK) == 0: |
---|
| 50 | mkdir (project.output_build_time_dir) |
---|
[3827] | 51 | copy (dirname(project.__file__) +sep+ project.__name__+'.py', |
---|
[3851] | 52 | project.output_build_time_dir + project.__name__+'.py') #copies project.py |
---|
| 53 | copy (__file__, project.output_build_time_dir + basename(__file__)) |
---|
| 54 | print 'files '+ project.__name__+'.py and '+ basename(__file__)+' copied to '+ project.output_build_time_dir #copies this file |
---|
| 55 | #import sys; sys.exit() |
---|
[3827] | 56 | |
---|
| 57 | #------------------------------------------------------------------------------- |
---|
| 58 | # Preparation of topographic data |
---|
| 59 | # |
---|
| 60 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 61 | # Do for coarse and fine data |
---|
| 62 | # Fine pts file to be clipped to area of interest |
---|
| 63 | #------------------------------------------------------------------------------- |
---|
| 64 | |
---|
| 65 | # topography directory filenames |
---|
| 66 | onshore_dir_name = project.onshore_dir_name |
---|
| 67 | coast_dir_name = project.coast_dir_name |
---|
| 68 | islands_dir_name = project.islands_dir_name |
---|
| 69 | offshore_dir_name = project.offshore_dir_name |
---|
| 70 | offshore_dir_name1 = project.offshore_dir_name1 |
---|
| 71 | offshore_dir_name2 = project.offshore_dir_name2 |
---|
| 72 | offshore_dir_name3 = project.offshore_dir_name3 |
---|
| 73 | offshore_dir_name4 = project.offshore_dir_name4 |
---|
| 74 | offshore_dir_name5 = project.offshore_dir_name5 |
---|
| 75 | offshore_dir_name6 = project.offshore_dir_name6 |
---|
| 76 | offshore_dir_name7 = project.offshore_dir_name7 |
---|
| 77 | offshore_dir_name8 = project.offshore_dir_name8 |
---|
| 78 | offshore_dir_name9 = project.offshore_dir_name9 |
---|
| 79 | offshore_dir_name10 = project.offshore_dir_name10 |
---|
| 80 | offshore_dir_name11 = project.offshore_dir_name11 |
---|
| 81 | offshore_dir_name12 = project.offshore_dir_name12 |
---|
| 82 | offshore_dir_name13 = project.offshore_dir_name13 |
---|
| 83 | offshore_dir_name14 = project.offshore_dir_name14 |
---|
| 84 | |
---|
| 85 | # creates DEM from asc data |
---|
| 86 | convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True) |
---|
| 87 | convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True) |
---|
| 88 | |
---|
| 89 | #creates pts file for onshore DEM |
---|
| 90 | dem2pts(onshore_dir_name, |
---|
| 91 | # easting_min=project.eastingmin, |
---|
| 92 | # easting_max=project.eastingmax, |
---|
| 93 | # northing_min=project.northingmin, |
---|
| 94 | # northing_max= project.northingmax, |
---|
| 95 | use_cache=True, |
---|
| 96 | verbose=True) |
---|
| 97 | |
---|
| 98 | #creates pts file for islands DEM |
---|
| 99 | dem2pts(islands_dir_name, use_cache=True, verbose=True) |
---|
| 100 | |
---|
[3851] | 101 | print'create Geospatial data objects from topographies' |
---|
[3827] | 102 | G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts') |
---|
| 103 | G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya') |
---|
| 104 | G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts') |
---|
| 105 | G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya') |
---|
| 106 | G_off1 = Geospatial_data(file_name = project.offshore_dir_name1 + '.xya') |
---|
| 107 | G_off2 = Geospatial_data(file_name = project.offshore_dir_name2 + '.xya') |
---|
| 108 | G_off3 = Geospatial_data(file_name = project.offshore_dir_name3 + '.xya') |
---|
| 109 | G_off4 = Geospatial_data(file_name = project.offshore_dir_name4 + '.xya') |
---|
| 110 | G_off5 = Geospatial_data(file_name = project.offshore_dir_name5 + '.xya') |
---|
| 111 | G_off6 = Geospatial_data(file_name = project.offshore_dir_name6 + '.xya') |
---|
| 112 | G_off7 = Geospatial_data(file_name = project.offshore_dir_name7 + '.xya') |
---|
| 113 | G_off8 = Geospatial_data(file_name = project.offshore_dir_name8 + '.xya') |
---|
| 114 | G_off9 = Geospatial_data(file_name = project.offshore_dir_name9 + '.xya') |
---|
| 115 | G_off10 = Geospatial_data(file_name = project.offshore_dir_name10 + '.xya') |
---|
| 116 | G_off11 = Geospatial_data(file_name = project.offshore_dir_name11 + '.xya') |
---|
| 117 | G_off12 = Geospatial_data(file_name = project.offshore_dir_name12 + '.xya') |
---|
| 118 | G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya') |
---|
| 119 | G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya') |
---|
| 120 | |
---|
[3851] | 121 | print'add all geospatial objects' |
---|
[3827] | 122 | G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \ |
---|
| 123 | + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \ |
---|
| 124 | + G_off13 + G_off14 |
---|
| 125 | |
---|
[3851] | 126 | print'clip combined geospatial object by bounding polygon' |
---|
[3827] | 127 | G.clip(project.bounding_polygon) |
---|
| 128 | #FIXME: add a clip function to pts |
---|
| 129 | |
---|
[3851] | 130 | print'export combined DEM file' |
---|
| 131 | if access(project.topographies_time_dir,F_OK) == 0: |
---|
| 132 | mkdir (project.topographies_time_dir) |
---|
| 133 | G.export_points_file(project.combined_time_dir_name + '.pts') |
---|
[3827] | 134 | |
---|
| 135 | #------------------------------------------------------------------------- |
---|
| 136 | # Convert URS to SWW file for boundary conditions |
---|
| 137 | #------------------------------------------------------------------------- |
---|
[3851] | 138 | print 'starting to create boundary conditions' |
---|
| 139 | boundaries_dir_name = project.boundaries_dir_name |
---|
[3827] | 140 | |
---|
[3851] | 141 | from anuga.shallow_water.data_manager import urs2sww |
---|
[3827] | 142 | |
---|
[3851] | 143 | print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary |
---|
| 144 | print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary |
---|
[3827] | 145 | |
---|
[3851] | 146 | if access(project.boundaries_time_dir,F_OK) == 0: |
---|
| 147 | mkdir (project.boundaries_time_dir) |
---|
| 148 | urs2sww(boundaries_dir_name,basename_out= project.boundaries_time_dir_name, |
---|
| 149 | minlat=project.south_boundary, maxlat=project.north_boundary, |
---|
[3827] | 150 | minlon= project.west_boundary, maxlon=project.east_boundary, |
---|
[3851] | 151 | mint=0, maxt= 35000, |
---|
[3827] | 152 | verbose='true') |
---|
[3851] | 153 | |
---|
[3827] | 154 | |
---|
| 155 | |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | |
---|
| 159 | |
---|
[3851] | 160 | |
---|
| 161 | |
---|