source: anuga_work/production/perth_2006/build_perth.py @ 4147

Last change on this file since 4147 was 4147, checked in by sexton, 18 years ago

(1) updates to Dampier script based on Perth script (2) minor updates to Onslow report

File size: 7.6 KB
Line 
1"""Script for running tsunami inundation scenario for Perth, WA, Australia.
2
3Source data such as elevation and boundary data is assumed to be available in
4directories specified by project.py
5The output sww file is stored in project.output_time_dir
6
7The scenario is defined by a triangular mesh created from project.polygon,
8the elevation data and a simulated submarine landslide.
9
10Ole 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
18from os import sep
19from os.path import dirname, basename
20from os import mkdir, access, F_OK
21from shutil import copy
22import time
23import sys
24
25# Related major packages
26from 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
30from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
31#from anuga.pmesh.mesh_interface import create_mesh_from_regions
32from anuga.geospatial_data.geospatial_data import *
33from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
34
35# Application specific imports
36import 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
43copy_code_files(project.output_build_time_dir,__file__, 
44               dirname(project.__file__)+sep+ project.__name__+'.py' )
45
46start_screen_catcher(project.output_build_time_dir)
47
48print '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#-------------------------------------------------------------------------------
57print"project.bounding_polygon",project.bounding_polygon
58print"project.combined_dir_name",project.combined_dir_name
59
60# topography directory filenames
61onshore_in_dir_name = project.onshore_in_dir_name
62coast_in_dir_name = project.coast_in_dir_name
63island_in_dir_name = project.island_in_dir_name
64island_in_dir_name1 = project.island_in_dir_name1
65island_in_dir_name2 = project.island_in_dir_name2
66island_in_dir_name3 = project.island_in_dir_name3
67offshore_in_dir_name = project.offshore_in_dir_name
68offshore1_in_dir_name = project.offshore1_in_dir_name
69
70onshore_dir_name = project.onshore_dir_name
71coast_dir_name = project.coast_dir_name
72island_dir_name = project.island_dir_name
73island_dir_name1 = project.island_dir_name1
74island_dir_name2 = project.island_dir_name2
75island_dir_name3 = project.island_dir_name3
76offshore_dir_name = project.offshore_dir_name
77
78# creates DEM from asc data
79print "creates DEMs from asc data"
80convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
81convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
82convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True)
83convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True)
84convert_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
87print "creates pts file for onshore DEM"
88dem2pts(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
97dem2pts(island_dir_name, use_cache=True, verbose=True)
98dem2pts(island_dir_name1, use_cache=True, verbose=True)
99dem2pts(island_dir_name2, use_cache=True, verbose=True)
100dem2pts(island_dir_name3, use_cache=True, verbose=True)
101
102print'create Geospatial data1 objects from topographies'
103G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
104print'create Geospatial data2 objects from topographies'
105G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
106print'create Geospatial data3 objects from topographies'
107G3 = Geospatial_data(file_name = island_dir_name + '.pts')
108print'create Geospatial data4 objects from topographies'
109G4 = Geospatial_data(file_name = island_dir_name1 + '.pts')
110print'create Geospatial data5 objects from topographies'
111G5 = Geospatial_data(file_name = island_dir_name2 + '.pts')
112print'create Geospatial data6 objects from topographies'
113G6 = Geospatial_data(file_name = island_dir_name3 + '.pts')
114print'create Geospatial data7 objects from topographies'
115G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
116print'create Geospatial data8 objects from topographies'
117G_off1 = Geospatial_data(file_name = offshore1_in_dir_name + '.xya')
118
119print'add all geospatial objects'
120G = G1 + G2 + G3 + G4 + G5 + G6 + G_off + G_off1
121
122print'clip combined geospatial object by bounding polygon'
123G_clipped = G.clip(project.bounding_polygon)
124#FIXME: add a clip function to pts
125#print'shape of clipped data', G_clipped.get_data_points().shape
126
127print'export combined DEM file'
128if access(project.topographies_dir,F_OK) == 0:
129    mkdir (project.topographies_dir)
130G_clipped.export_points_file(project.combined_dir_name + '.pts')
131#G_clipped.export_points_file(project.combined_dir_name + '.xya')
132
133'''
134print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
135G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
136print'split'
137G_all_1, G_all_2 = G_all.split(.10)
138print'export 1'
139G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
140print'export 2'
141G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
142
143
144#-------------------------------------------------------------------------
145# Convert URS to SWW file for boundary conditions
146#-------------------------------------------------------------------------
147print 'starting to create boundary conditions'
148boundaries_in_dir_name = project.boundaries_in_dir_name
149
150from anuga.shallow_water.data_manager import urs2sww
151
152print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
153print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
154
155#import sys; sys.exit()
156
157#if access(project.boundaries_dir,F_OK) == 0:
158#    mkdir (project.boundaries_dir)
159
160from caching import cache
161cache(urs2sww,
162      (boundaries_in_dir_name,
163       project.boundaries_dir_name1),
164      {'verbose': True,
165       'minlat': project.south_boundary,
166       'maxlat': project.north_boundary,
167       'minlon': project.west_boundary,
168       'maxlon': project.east_boundary,
169#       'minlat': project.south,
170#       'maxlat': project.north,
171#       'minlon': project.west,
172#       'maxlon': project.east,
173       'mint': 0, 'maxt': 40000,
174#       'origin': domain.geo_reference.get_origin(),
175       'mean_stage': project.tide,
176#       'zscale': 1,                 #Enhance tsunami
177       'fail_on_NaN': False},
178       verbose = True,
179       )
180#       dependencies = source_dir + project.boundary_basename + '.sww')
181
182'''
183
184
185
186
187
188
189
Note: See TracBrowser for help on using the repository browser.