source: anuga_work/production/busselton/build_busselton.py @ 5381

Last change on this file since 5381 was 5381, checked in by kristy, 16 years ago

update scripts

File size: 5.9 KB
RevLine 
[4132]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
28from anuga.geospatial_data.geospatial_data import *
[5149]29from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
[4132]30
31# Application specific imports
32import project   # Definition of file names and polygons
33
34#------------------------------------------------------------------------------
35# Copy scripts to time stamped output directory and capture screen
36# output to file
37#------------------------------------------------------------------------------
38
39start_screen_catcher(project.output_build_time_dir)
40
[5149]41print 'time stamp: ',project.time
[4132]42print 'USER: ', project.user
43
44
45copy_code_files(project.output_build_time_dir,__file__, 
46               dirname(project.__file__)+sep+ project.__name__+'.py' )
47
48
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#-------------------------------------------------------------------------------
57
58# topography directory filenames
59onshore_in_dir_name = project.onshore_in_dir_name
60coast_in_dir_name = project.coast_in_dir_name
[5149]61coast_in_dir_name1 = project.coast_in_dir_name1
[4132]62#island_in_dir_name = project.island_in_dir_name
63offshore_in_dir_name = project.offshore_in_dir_name
64
65onshore_dir_name = project.onshore_dir_name
66coast_dir_name = project.coast_dir_name
[5149]67coast_dir_name1 = project.coast_dir_name1
[4132]68#island_dir_name = project.island_dir_name
69offshore_dir_name = project.offshore_dir_name
70
71# creates DEM from asc data
72print "creates DEMs from ascii data"
73convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
74
75#creates pts file for onshore DEM
76print "creates pts file for onshore DEM"
77dem2pts(onshore_dir_name,
78#        easting_min=project.eastingmin,
79#        easting_max=project.eastingmax,
80#        northing_min=project.northingmin,
81#        northing_max= project.northingmax,
82        use_cache=True, 
83        verbose=True)
84
85#creates pts file for island DEM
86#dem2pts(island_dir_name, use_cache=True, verbose=True)
87
88print'create Geospatial data1 objects from topographies'
[5149]89G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
90G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
91G3 = Geospatial_data(file_name = coast_in_dir_name1 + '.txt',verbose=True)
92#G3 = Geospatial_data(file_name = island_dir_name + '.pts',verbose=True)
93G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
[4132]94
95print'add all geospatial objects'
[5149]96G = G1 + G2 + G3 + G_off
[4132]97
98print'clip combined geospatial object by bounding polygon'
99#G_clipped = G.clip(project.bounding_polygon)
100#FIXME: add a clip function to pts
101#print'shape of clipped data', G_clipped.get_data_points().shape
102
103print'export combined DEM file'
104if access(project.topographies_dir,F_OK) == 0:
105    mkdir (project.topographies_dir)
[5149]106G.export_points_file(project.combined_dir_name + '.txt')
[4132]107#G_clipped.export_points_file(project.combined_dir_name + '.xya')
108
[5149]109print'project.combined_dir_name + .txt',project.combined_dir_name + '.txt'
110#G_all=Geospatial_data(file_name = project.combined_dir_name + '.txt')
[4132]111print'split'
[5381]112G_all_1 = G.split(.10)
[4132]113print'export 1'
[5149]114G_all_1.export_points_file(project.combined_dir_name+'_small' + '.txt')
115#print'export 2'
116#G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
[4132]117
[5149]118'''
[4132]119#-------------------------------------------------------------------------
120# Convert URS to SWW file for boundary conditions
121#-------------------------------------------------------------------------
122print 'starting to create boundary conditions'
123boundaries_in_dir_name = project.boundaries_in_dir_name
124
125from anuga.shallow_water.data_manager import urs2sww
126
127print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
128print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
129
130#import sys; sys.exit()
131
132#if access(project.boundaries_dir,F_OK) == 0:
133#    mkdir (project.boundaries_dir)
134
135from caching import cache
136cache(urs2sww,
137      (boundaries_in_dir_name,
138       project.boundaries_dir_name1),
139      {'verbose': True,
140       'minlat': project.south_boundary,
141       'maxlat': project.north_boundary,
142       'minlon': project.west_boundary,
143       'maxlon': project.east_boundary,
144#       'minlat': project.south,
145#       'maxlat': project.north,
146#       'minlon': project.west,
147#       'maxlon': project.east,
148       'mint': 0, 'maxt': 40000,
149#       'origin': domain.geo_reference.get_origin(),
150       'mean_stage': project.tide,
151#       'zscale': 1,                 #Enhance tsunami
152       'fail_on_NaN': False},
153       verbose = True,
154       )
155#       dependencies = source_dir + project.boundary_basename + '.sww')
156
157'''
158
159
160
161
162
163
164
Note: See TracBrowser for help on using the repository browser.