source: anuga_work/production/dampier_2006/build_dampier.py @ 4430

Last change on this file since 4430 was 4357, checked in by nick, 17 years ago

update to dampier

File size: 6.9 KB
Line 
1"""Script for running tsunami inundation scenario for Dampier, 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
26# Related major packages
27from anuga.shallow_water import Domain
28from anuga.shallow_water import Dirichlet_boundary
29from anuga.shallow_water import File_boundary
30from anuga.shallow_water import Reflective_boundary
31from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
32from anuga.pmesh.mesh_interface import create_mesh_from_regions
33from anuga.geospatial_data.geospatial_data import *
34from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
35from anuga_parallel.parallel_abstraction import get_processor_name
36
37# Application specific imports
38import project   # Definition of file names and polygons
39
40#------------------------------------------------------------------------------
41# Copy scripts to time stamped output directory and capture screen
42# output to file
43#------------------------------------------------------------------------------
44copy_code_files(project.output_build_time_dir,__file__, 
45               dirname(project.__file__)+sep+ project.__name__+'.py' )
46
47start_screen_catcher(project.output_build_time_dir)
48print "Processor Name:",get_processor_name()
49
50print 'USER:    ', project.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#-------------------------------------------------------------------------------
59print"project.combined_dir_name",project.combined_dir_name
60
61# topography directory filenames
62onshore_in_dir_name = project.onshore_in_dir_name
63coast_in_dir_name = project.coast_in_dir_name
64offshore_in_dir_name = project.offshore_in_dir_name
65offshore1_in_dir_name = project.offshore1_in_dir_name
66offshore2_in_dir_name = project.offshore2_in_dir_name
67
68onshore_dir_name = project.onshore_dir_name
69coast_dir_name = project.coast_dir_name
70offshore_dir_name = project.offshore_dir_name
71offshore1_dir_name = project.offshore1_dir_name
72offshore2_dir_name = project.offshore2_dir_name
73'''
74# creates DEM from asc data
75print "creates DEMs from asc data"
76convert_dem_from_ascii2netcdf(onshore_in_dir_name,
77                              basename_out=onshore_dir_name,
78                              use_cache=True, verbose=True)
79
80convert_dem_from_ascii2netcdf(offshore1_in_dir_name,
81                              basename_out=offshore1_dir_name,
82                              use_cache=True, verbose=True)
83convert_dem_from_ascii2netcdf(offshore2_in_dir_name,
84                              basename_out=offshore2_dir_name,
85                              use_cache=True, verbose=True)
86
87#creates pts file for onshore DEM
88print "creates pts file for onshore DEM"
89dem2pts(onshore_dir_name, use_cache=True, verbose=True)
90
91#creates pts file for clipped DEMs
92dem2pts(offshore1_dir_name,
93        basename_out=offshore1_dir_name,
94        use_cache=True, verbose=True)
95dem2pts(offshore2_dir_name,
96        basename_out=offshore2_dir_name,
97        use_cache=True, verbose=True)
98
99print'create Geospatial onshore objects from topographies'
100#print'create Geospatial coastal objects from topographies'
101#print'create Geospatial nickel bay objects from topographies'
102#print'create Geospatial clipped TIN offshore objects from topographies'
103#print'create Geospatial offshore objects from topographies'
104G = Geospatial_data(file_name = onshore_dir_name + '.pts') +\
105    Geospatial_data(file_name = coast_in_dir_name + '.txt') +\
106    Geospatial_data(file_name = offshore_in_dir_name + '.txt') +\
107    Geospatial_data(file_name = offshore1_dir_name + '.pts') +\
108    Geospatial_data(file_name = offshore2_dir_name + '.pts')
109     
110
111print'add all geospatial objects'
112#G = G1 + G2 + G3 + G4 + G_off
113
114print'clip combined geospatial object by bounding polygon'
115G_clipped = G.clip(project.poly_all)
116#FIXME: add a clip function to pts
117#print'shape of clipped data', G_clipped.get_data_points().shape
118
119print'export combined DEM file'
120if access(project.topographies_dir,F_OK) == 0:
121    mkdir (project.topographies_dir)
122G_clipped.export_points_file(project.combined_dir_name + '.pts')
123#G_clipped.export_points_file(project.combined_dir_name + '.txt')
124
125print'split combined data set'
126G_small, G_other = G_clipped.split(0.1, True)
127
128print'export split DEM file'
129G_small.export_points_file(project.combined_dir_name + '_small' + '.pts')
130G_other.export_points_file(project.combined_dir_name + '_other' + '.pts')
131
132print 'start reading:',project.combined_smaller_dir_name + '.txt'
133G = Geospatial_data(file_name = project.combined_smaller_dir_name + '.txt')
134print 'start split'
135G_smallest, G_other = G.split(0.1,True)
136
137print 'start export',project.combined_smallest_dir_name + '.txt'
138#G.export_points_file(project.combined_smaller_dir_name + '.txt')
139#G_smallest.export_points_file(project.combined_smallest_dir_name + '.txt')
140'''
141#-------------------------------------------------------------------------
142# Convert URS to SWW file for boundary conditions
143#-------------------------------------------------------------------------
144print 'starting to create boundary conditions'
145boundaries_in_dir_name = project.boundaries_in_dir_name
146
147from anuga.shallow_water.data_manager import urs2sww
148
149print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
150print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
151
152#import sys; sys.exit()
153
154#if access(project.boundaries_dir,F_OK) == 0:
155#    mkdir (project.boundaries_dir)
156
157from caching import cache
158cache(urs2sww,
159      (boundaries_in_dir_name,
160       project.boundaries_dir_name+'a'),
161      {'verbose': True,
162       'minlat': project.south_boundary,
163       'maxlat': project.north_boundary,
164       'minlon': project.west_boundary,
165       'maxlon': project.east_boundary,
166       'mint': 5000, 'maxt': 35000,
167#       'origin': domain.geo_reference.get_origin(),
168#       'mean_stage': project.tide,
169       'zscale': 1,                 #Enhance tsunami
170       'fail_on_NaN': False},
171       verbose = True,
172       )
173#       dependencies = source_dir + project.boundary_basename + '.sww')
174
175
176
177
178
Note: See TracBrowser for help on using the repository browser.