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_urs.py |
---|
5 | The output sww file is stored in project_urs.output_time_dir |
---|
6 | |
---|
7 | The scenario is defined by a triangular mesh created from project_urs.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 Geospatial_data |
---|
33 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files |
---|
34 | from anuga_parallel.parallel_abstraction import get_processor_name |
---|
35 | |
---|
36 | # Application specific imports |
---|
37 | import project_urs # Definition of file names and polygons |
---|
38 | |
---|
39 | #------------------------------------------------------------------------------ |
---|
40 | # Copy scripts to time stamped output directory and capture screen |
---|
41 | # output to file |
---|
42 | #------------------------------------------------------------------------------ |
---|
43 | |
---|
44 | copy_code_files(project_urs.output_build_time_dir,__file__, |
---|
45 | dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' ) |
---|
46 | |
---|
47 | start_screen_catcher(project_urs.output_build_time_dir) |
---|
48 | |
---|
49 | print "Processor Name:",get_processor_name() |
---|
50 | print 'USER: ', project_urs.user |
---|
51 | |
---|
52 | #------------------------------------------------------------------------------- |
---|
53 | # Preparation of topographic data |
---|
54 | # |
---|
55 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
56 | # Do for coarse and fine data |
---|
57 | # Fine pts file to be clipped to area of interest |
---|
58 | #------------------------------------------------------------------------------- |
---|
59 | print"project_urs.combined_dir_name",project_urs.combined_dir_name |
---|
60 | ''' |
---|
61 | # topography directory filenames |
---|
62 | onshore_in_dir_name = project_urs.onshore_in_dir_name |
---|
63 | coast_in_dir_name = project_urs.coast_in_dir_name |
---|
64 | offshore_in_dir_name = project_urs.offshore_in_dir_name |
---|
65 | offshore1_in_dir_name = project_urs.offshore1_in_dir_name |
---|
66 | offshore2_in_dir_name = project_urs.offshore2_in_dir_name |
---|
67 | |
---|
68 | onshore_dir_name = project_urs.onshore_dir_name |
---|
69 | coast_dir_name = project_urs.coast_dir_name |
---|
70 | offshore_dir_name = project_urs.offshore_dir_name |
---|
71 | offshore1_dir_name = project_urs.offshore1_dir_name |
---|
72 | offshore2_dir_name = project_urs.offshore2_dir_name |
---|
73 | |
---|
74 | # creates DEM from asc data |
---|
75 | print "creates DEMs from asc data" |
---|
76 | convert_dem_from_ascii2netcdf(onshore_in_dir_name, |
---|
77 | basename_out=onshore_dir_name, |
---|
78 | use_cache=True, verbose=True) |
---|
79 | #creates pts file for onshore DEM |
---|
80 | print "creates pts file for onshore DEM" |
---|
81 | dem2pts(onshore_dir_name, use_cache=True, verbose=True) |
---|
82 | |
---|
83 | print'create Geospatial onshore objects from topographies' |
---|
84 | G = Geospatial_data(file_name = onshore_dir_name + '.pts') +\ |
---|
85 | Geospatial_data(file_name = coast_in_dir_name + '.txt') +\ |
---|
86 | Geospatial_data(file_name = offshore_in_dir_name + '.txt') |
---|
87 | # +\ |
---|
88 | # Geospatial_data(file_name = offshore1_dir_name + '.pts') +\ |
---|
89 | # Geospatial_data(file_name = offshore2_dir_name + '.pts') |
---|
90 | |
---|
91 | |
---|
92 | print'clip combined geospatial object by bounding polygon' |
---|
93 | G_clipped = G.clip(project_urs.poly_all) |
---|
94 | |
---|
95 | print'export combined DEM file' |
---|
96 | if access(project_urs.topographies_dir,F_OK) == 0: |
---|
97 | mkdir (project_urs.topographies_dir) |
---|
98 | G.export_points_file(project_urs.combined_dir_name + '.txt') |
---|
99 | |
---|
100 | print'split combined data set' |
---|
101 | G_small, G_other = G_clipped.split(0.1, True) |
---|
102 | |
---|
103 | print'export split DEM file' |
---|
104 | G_small.export_points_file(project_urs.combined_dir_name + '_small' + '.txt') |
---|
105 | #G_other.export_points_file(project_urs.combined_dir_name + '_other' + '.pts') |
---|
106 | |
---|
107 | ''' |
---|
108 | print 'start reading:',project_urs.combined_dir_name + '.txt' |
---|
109 | G = Geospatial_data(file_name = (project_urs.combined_dir_name + '.txt')) |
---|
110 | G_clipped = G.clip(project_urs.poly_topo_clip) |
---|
111 | |
---|
112 | print 'start split' |
---|
113 | #G_smallest, G_other = G.split(0.1,True) |
---|
114 | |
---|
115 | print 'start export',project_urs.combined_smallest_dir_name + '.txt' |
---|
116 | #G.export_points_file(project_urs.combined_smaller_dir_name + '.txt') |
---|
117 | G_clipped.export_points_file(project_urs.combined_smallest_dir_name + '.txt') |
---|
118 | |
---|
119 | #------------------------------------------------------------------------- |
---|
120 | # Convert URS to SWW file for boundary conditions |
---|
121 | #------------------------------------------------------------------------- |
---|
122 | print 'starting to create boundary conditions' |
---|
123 | #boundaries_in_dir_name = project_urs.boundaries_in_dir_name |
---|
124 | |
---|
125 | #print 'minlat=project_urs.north_boundary, maxlat=project_urs.south_boundary',project_urs.north_boundary, project_urs.south_boundary |
---|
126 | #print 'minlon= project_urs.west_boundary, maxlon=project_urs.east_boundary',project_urs.west_boundary, project_urs.east_boundary |
---|
127 | ''' |
---|
128 | from anuga.shallow_water.data_manager import urs_ungridded2sww |
---|
129 | |
---|
130 | print 'boundaries_in_dir_name',project_urs.boundaries_in_dir_name |
---|
131 | print 'project.boundaries_dir_name',project_urs.boundaries_dir_name |
---|
132 | |
---|
133 | urs_ungridded2sww(project_urs.boundaries_in_dir_name, project_urs.boundaries_dir_name, |
---|
134 | verbose=True, mint=5000, maxt=35000, zscale=1) |
---|
135 | ''' |
---|
136 | |
---|
137 | |
---|
138 | |
---|
139 | |
---|