1 | """Script for running tsunami inundation scenario for Perth, 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 |
---|
27 | #from anuga.shallow_water import Dirichlet_boundary |
---|
28 | #from anuga.shallow_water import File_boundary |
---|
29 | #from anuga.shallow_water import Reflective_boundary |
---|
30 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
31 | #from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
32 | from anuga.geospatial_data.geospatial_data import * |
---|
33 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
---|
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 | #------------------------------------------------------------------------------ |
---|
42 | |
---|
43 | copy_code_files(project.output_build_time_dir,__file__, |
---|
44 | dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
45 | |
---|
46 | start_screen_catcher(project.output_build_time_dir) |
---|
47 | |
---|
48 | print 'USER: ', project.user |
---|
49 | |
---|
50 | #------------------------------------------------------------------------------- |
---|
51 | # Preparation of topographic data |
---|
52 | # |
---|
53 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
54 | # Do for coarse and fine data |
---|
55 | # Fine pts file to be clipped to area of interest |
---|
56 | #------------------------------------------------------------------------------- |
---|
57 | print"project.poly_all",project.poly_all |
---|
58 | print"project.combined_dir_name",project.combined_dir_name |
---|
59 | |
---|
60 | # topography directory filenames |
---|
61 | onshore_in_dir_name = project.onshore_in_dir_name |
---|
62 | coast_in_dir_name = project.coast_in_dir_name |
---|
63 | island_in_dir_name = project.island_in_dir_name |
---|
64 | island_in_dir_name1 = project.island_in_dir_name1 |
---|
65 | island_in_dir_name2 = project.island_in_dir_name2 |
---|
66 | island_in_dir_name3 = project.island_in_dir_name3 |
---|
67 | offshore_in_dir_name = project.offshore_in_dir_name |
---|
68 | |
---|
69 | |
---|
70 | onshore_dir_name = project.onshore_dir_name |
---|
71 | coast_dir_name = project.coast_dir_name |
---|
72 | island_dir_name = project.island_dir_name |
---|
73 | island_dir_name1 = project.island_dir_name1 |
---|
74 | island_dir_name2 = project.island_dir_name2 |
---|
75 | island_dir_name3 = project.island_dir_name3 |
---|
76 | offshore_dir_name = project.offshore_dir_name |
---|
77 | |
---|
78 | # creates DEM from asc data |
---|
79 | print "creates DEMs from asc data" |
---|
80 | convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True) |
---|
81 | convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True) |
---|
82 | convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True) |
---|
83 | convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True) |
---|
84 | convert_dem_from_ascii2netcdf(island_in_dir_name3, basename_out=island_dir_name3, use_cache=True, verbose=True) |
---|
85 | |
---|
86 | #creates pts file for onshore DEM |
---|
87 | print "creates pts file for onshore DEM" |
---|
88 | dem2pts(onshore_dir_name, |
---|
89 | # easting_min=project.eastingmin, |
---|
90 | # easting_max=project.eastingmax, |
---|
91 | # northing_min=project.northingmin, |
---|
92 | # northing_max= project.northingmax, |
---|
93 | use_cache=True, |
---|
94 | verbose=True) |
---|
95 | |
---|
96 | #creates pts file for island DEM |
---|
97 | dem2pts(island_dir_name, use_cache=True, verbose=True) |
---|
98 | dem2pts(island_dir_name1, use_cache=True, verbose=True) |
---|
99 | dem2pts(island_dir_name2, use_cache=True, verbose=True) |
---|
100 | dem2pts(island_dir_name3, use_cache=True, verbose=True) |
---|
101 | |
---|
102 | print'create Geospatial data1 objects from topographies' |
---|
103 | G1 = Geospatial_data(file_name = onshore_dir_name + '.pts') |
---|
104 | print'create Geospatial data2 objects from topographies' |
---|
105 | G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt') |
---|
106 | print'create Geospatial data3 objects from topographies' |
---|
107 | G3 = Geospatial_data(file_name = island_dir_name + '.pts') |
---|
108 | print'create Geospatial data4 objects from topographies' |
---|
109 | G4 = Geospatial_data(file_name = island_dir_name1 + '.pts') |
---|
110 | print'create Geospatial data5 objects from topographies' |
---|
111 | G5 = Geospatial_data(file_name = island_dir_name2 + '.pts') |
---|
112 | print'create Geospatial data6 objects from topographies' |
---|
113 | G6 = Geospatial_data(file_name = island_dir_name3 + '.pts') |
---|
114 | print'create Geospatial data7 objects from topographies' |
---|
115 | G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt') |
---|
116 | |
---|
117 | |
---|
118 | print'add all geospatial objects' |
---|
119 | G = G1 + G2 + G3 + G4 + G5 + G6 + G_off |
---|
120 | |
---|
121 | print'clip combined geospatial object by bounding polygon' |
---|
122 | G_clipped = G.clip(project.poly_all) |
---|
123 | |
---|
124 | print'export combined DEM file' |
---|
125 | if access(project.topographies_dir,F_OK) == 0: |
---|
126 | mkdir (project.topographies_dir) |
---|
127 | G_clipped.export_points_file(project.combined_dir_name + '.pts') |
---|
128 | #G_clipped.export_points_file(project.combined_dir_name + '.xya') - use for comparision in ARC |
---|
129 | |
---|
130 | import sys |
---|
131 | sys.exit() |
---|
132 | |
---|
133 | #------------------------------------------------------------------------- |
---|
134 | # Convert URS to SWW file for boundary conditions |
---|
135 | #------------------------------------------------------------------------- |
---|
136 | print 'starting to create boundary conditions' |
---|
137 | from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww |
---|
138 | |
---|
139 | boundaries_in_dir_name = project.boundaries_in_dir_name |
---|
140 | print 'boundaries_in_dir_name',project.boundaries_in_dir_name |
---|
141 | |
---|
142 | |
---|
143 | #import sys; sys.exit() |
---|
144 | |
---|
145 | urs_ungridded2sww(project.boundaries_in_dir_name, project.boundaries_in_dir_name, |
---|
146 | verbose=True, mint=4000, maxt=80000, zscale=1) |
---|
147 | |
---|
148 | |
---|
149 | |
---|
150 | |
---|
151 | |
---|
152 | |
---|
153 | |
---|
154 | |
---|
155 | |
---|
156 | |
---|