source: anuga_work/production/geraldton/build_geraldton.py @ 5150

Last change on this file since 5150 was 5150, checked in by nick, 16 years ago

update geraldton

File size: 7.0 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# Related major packages
26from anuga.shallow_water import Domain
27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,start_screen_catcher, copy_code_files
28from anuga.geospatial_data.geospatial_data import *
29
30# Application specific imports
31import project   # Definition of file names and polygons
32
33#------------------------------------------------------------------------------
34# Copy scripts to time stamped output directory and capture screen
35# output to file
36#------------------------------------------------------------------------------
37
38print 'time stamp: ',project.build_time
39
40copy_code_files(project.output_build_time_dir,__file__, 
41               dirname(project.__file__)+sep+ project.__name__+'.py' )
42
43start_screen_catcher(project.output_build_time_dir)
44
45print 'time stamp: ',project.build_time
46print 'USER: ', project.user, project.host
47
48
49#-------------------------------------------------------------------------------
50# Preparation of topographic data
51#
52# Convert ASC 2 DEM 2 PTS using source data and store result in source data
53# Do for coarse and fine data
54# Fine pts file to be clipped to area of interest
55#-------------------------------------------------------------------------------
56
57# topography directory filenames
58onshore_in_dir_name = project.onshore_in_dir_name
59coast_in_dir_name = project.coast_in_dir_name
60island_in_dir_name = project.island_in_dir_name
61offshore_in_dir_name = project.offshore_in_dir_name
62
63onshore_dir_name = project.onshore_dir_name
64coast_dir_name = project.coast_dir_name
65island_dir_name = project.island_dir_name
66offshore_dir_name = project.offshore_dir_name
67
68# creates DEM from asc data
69print "creates DEMs from ascii data"
70convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
71convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
72
73#creates pts file for onshore DEM
74print "creates pts file for onshore DEM"
75dem2pts(onshore_dir_name,
76        use_cache=True, 
77        verbose=True)
78
79#creates pts file for island DEM
80dem2pts(island_dir_name, use_cache=True, verbose=True)
81
82print'create Geospatial data1 objects from topographies'
83G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
84print'finished reading in %s.pts' %onshore_dir_name
85G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
86G3 = Geospatial_data(file_name = island_dir_name + '.pts',verbose=True)
87G4 = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
88print'finished reading in files'
89
90#G_onshore_clip = G1.clip(project.poly_all, verbose=True)
91#print 'finished clip'
92#G_onshore_clip_clipout=G_onshore_clip.clip_outside(project.poly_cbd, verbose=True)
93#print 'finished clipout'
94## reduce resolution by 10 time (1/10) of the original file. This was determined
95## as the resolution for this region is going to be 20000 and the asc grid cellsize is
96## 10m there (100). There is 2 20000m triangles in a 200m x 200m grid sqrt of 40000 is about 200,
97## a 200m x 200m area would contain 400 points at a 10m grid spacing and we only need 4. So 4/400 is 1/100 or 0.01
98## 1/100 reduction is possible but a 1/10 is enough to help with file size and loading
99#G_onshore_clip_clipout_reduced=G_onshore_clip_clipout.split(0.1, verbose=True)
100#
101#G_cbd = G1.clip(projec.poly_cbd,verbose=True)
102##G_cbd_reduced = G_cbd.split(0.25, verbose=True)
103#
104#G_bathy = G4.clip(project.poly_all)
105#G_island = G3.clip(project.poly_all)
106#G_coast = G2.clip(project.poly_all)
107#
108#print'add all geospatial objects'
109#G = G_onshore_clip_clipout_reduced + G_cbd_reduced + G_bathy + G_island + G_coast
110
111G = G1 + G2 + G3 + G4
112
113G_clip = G.clip(project.poly_all, verbose=True)
114
115print'clip combined geospatial object by bounding polygon'
116#G_clipped = G.clip(project.bounding_polygon)
117#FIXME: add a clip function to pts
118#print'shape of clipped data', G_clipped.get_data_points().shape
119
120print'export combined DEM file'
121if access(project.topographies_dir,F_OK) == 0:
122    mkdir (project.topographies_dir)
123G_clip.export_points_file(project.combined_dir_name + '.txt')
124#G_clipped.export_points_file(project.combined_dir_name + '.xya')
125G_small,G_other = G_clip.split(0.01, verbose=True)
126G_small.export_points_file(project.combined_dir_name_small + '.txt')
127
128'''
129print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
130G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
131print'split'
132G_all_1, G_all_2 = G_all.split(.10)
133print'export 1'
134G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
135print'export 2'
136G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
137
138
139#-------------------------------------------------------------------------
140# Convert URS to SWW file for boundary conditions
141#-------------------------------------------------------------------------
142print 'starting to create boundary conditions'
143boundaries_in_dir_name = project.boundaries_in_dir_name
144
145from anuga.shallow_water.data_manager import urs2sww
146
147print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
148print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
149
150#import sys; sys.exit()
151
152#if access(project.boundaries_dir,F_OK) == 0:
153#    mkdir (project.boundaries_dir)
154
155from caching import cache
156cache(urs2sww,
157      (boundaries_in_dir_name,
158       project.boundaries_dir_name1),
159      {'verbose': True,
160       'minlat': project.south_boundary,
161       'maxlat': project.north_boundary,
162       'minlon': project.west_boundary,
163       'maxlon': project.east_boundary,
164#       'minlat': project.south,
165#       'maxlat': project.north,
166#       'minlon': project.west,
167#       'maxlon': project.east,
168       'mint': 0, 'maxt': 40000,
169#       'origin': domain.geo_reference.get_origin(),
170       'mean_stage': project.tide,
171#       'zscale': 1,                 #Enhance tsunami
172       'fail_on_NaN': False},
173       verbose = True,
174       )
175#       dependencies = source_dir + project.boundary_basename + '.sww')
176
177'''
178
179
180
181
182
183
184
Note: See TracBrowser for help on using the repository browser.