source: anuga_work/production/perth/build_perth.py @ 5524

Last change on this file since 5524 was 5524, checked in by kristy, 14 years ago

input more bathymetry data

File size: 6.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.shallow_water.data_manager 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.poly_all",project.poly_all
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
68offshore_in_dir_name1 = project.offshore_in_dir_name1
69offshore_in_dir_name2 = project.offshore_in_dir_name2
70
71
72onshore_dir_name = project.onshore_dir_name
73coast_dir_name = project.coast_dir_name
74island_dir_name = project.island_dir_name
75island_dir_name1 = project.island_dir_name1
76island_dir_name2 = project.island_dir_name2
77island_dir_name3 = project.island_dir_name3
78offshore_dir_name = project.offshore_dir_name
79offshore_dir_name1 = project.offshore_dir_name1
80offshore_dir_name2 = project.offshore_dir_name2
81
82# creates DEM from asc data
83print "creates DEMs from asc data"
84convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
85convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
86convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True)
87convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True)
88convert_dem_from_ascii2netcdf(island_in_dir_name3, basename_out=island_dir_name3, use_cache=True, verbose=True)
89
90#creates pts file for onshore DEM
91print "creates pts file for onshore DEM"
92dem2pts(onshore_dir_name,
93#        easting_min=project.eastingmin,
94#        easting_max=project.eastingmax,
95#        northing_min=project.northingmin,
96#        northing_max= project.northingmax,
97        use_cache=True, 
98        verbose=True)
99
100#creates pts file for island DEM
101dem2pts(island_dir_name, use_cache=True, verbose=True)
102dem2pts(island_dir_name1, use_cache=True, verbose=True)
103dem2pts(island_dir_name2, use_cache=True, verbose=True)
104dem2pts(island_dir_name3, use_cache=True, verbose=True)
105
106print'create Geospatial data1 objects from topographies'
107G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
108print'create Geospatial data2 objects from topographies'
109G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt')
110print'create Geospatial data3 objects from topographies'
111G3 = Geospatial_data(file_name = island_dir_name + '.pts')
112print'create Geospatial data4 objects from topographies'
113G4 = Geospatial_data(file_name = island_dir_name1 + '.pts')
114print'create Geospatial data5 objects from topographies'
115G5 = Geospatial_data(file_name = island_dir_name2 + '.pts')
116print'create Geospatial data6 objects from topographies'
117G6 = Geospatial_data(file_name = island_dir_name3 + '.pts')
118print'create Geospatial data7 objects from topographies'
119G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt')
120print'create Geospatial data7 objects from topographies'
121G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt')
122print'create Geospatial data7 objects from topographies'
123G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt')
124
125
126print'add all geospatial objects'
127G = G1 + G2 + G3 + G4 + G5 + G6 + G_off + G_off1 + G_off2
128
129print'clip combined geospatial object by bounding polygon'
130G_clipped = G.clip(project.poly_all)
131
132print'export combined DEM file'
133if access(project.topographies_dir,F_OK) == 0:
134    mkdir (project.topographies_dir)
135G_clipped.export_points_file(project.combined_dir_name + '.pts')
136#G_clipped.export_points_file(project.combined_dir_name + '.xya') - use for comparision in ARC
137
138import sys
139sys.exit()
140
141#-------------------------------------------------------------------------
142# Convert URS to SWW file for boundary conditions
143#-------------------------------------------------------------------------
144print 'starting to create boundary conditions'
145from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
146
147boundaries_in_dir_name = project.boundaries_in_dir_name
148print 'boundaries_in_dir_name',project.boundaries_in_dir_name
149
150
151#import sys; sys.exit()
152
153urs_ungridded2sww(project.boundaries_in_dir_name, project.boundaries_in_dir_name,
154                  verbose=True, mint=4000, maxt=80000, zscale=1)
155
156
157
158
159
160
161
162
163
164
Note: See TracBrowser for help on using the repository browser.