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