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

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

update to dampier (this works)

File size: 6.9 KB
RevLine 
[3827]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 *
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#------------------------------------------------------------------------------
[3853]42
[3851]43if access(project.output_build_time_dir,F_OK) == 0:
44    mkdir (project.output_build_time_dir)
[3827]45copy (dirname(project.__file__) +sep+ project.__name__+'.py',
[3851]46      project.output_build_time_dir + project.__name__+'.py')     #copies project.py
47copy (__file__, project.output_build_time_dir + basename(__file__))
48print 'files '+ project.__name__+'.py and '+ basename(__file__)+' copied to '+ project.output_build_time_dir       #copies this file
49#import sys; sys.exit()
[3827]50
[3853]51
52#normal screen output is stored in
53screen_output_name = project.outputtimedir + "screen_output.txt"
54screen_error_name = project.outputtimedir + "screen_error.txt"
55
56#used to catch screen output to file
57sys.stdout = Screen_Catcher(screen_output_name)
58sys.stderr = Screen_Catcher(screen_error_name)
59
60print 'USER:    ', project.user
61
[3827]62#-------------------------------------------------------------------------------
63# Preparation of topographic data
64#
65# Convert ASC 2 DEM 2 PTS using source data and store result in source data
66# Do for coarse and fine data
67# Fine pts file to be clipped to area of interest
68#-------------------------------------------------------------------------------
69
70# topography directory filenames
71onshore_dir_name = project.onshore_dir_name
72coast_dir_name = project.coast_dir_name
73islands_dir_name = project.islands_dir_name
74offshore_dir_name = project.offshore_dir_name
75offshore_dir_name1 = project.offshore_dir_name1
76offshore_dir_name2 = project.offshore_dir_name2
77offshore_dir_name3 = project.offshore_dir_name3
78offshore_dir_name4 = project.offshore_dir_name4
79offshore_dir_name5 = project.offshore_dir_name5
80offshore_dir_name6 = project.offshore_dir_name6
81offshore_dir_name7 = project.offshore_dir_name7
82offshore_dir_name8 = project.offshore_dir_name8
83offshore_dir_name9 = project.offshore_dir_name9
84offshore_dir_name10 = project.offshore_dir_name10
85offshore_dir_name11 = project.offshore_dir_name11
86offshore_dir_name12 = project.offshore_dir_name12
87offshore_dir_name13 = project.offshore_dir_name13
88offshore_dir_name14 = project.offshore_dir_name14
89
90# creates DEM from asc data
91convert_dem_from_ascii2netcdf(onshore_dir_name, use_cache=True, verbose=True)
92convert_dem_from_ascii2netcdf(islands_dir_name, use_cache=True, verbose=True)
93
94#creates pts file for onshore DEM
95dem2pts(onshore_dir_name,
96#        easting_min=project.eastingmin,
97#        easting_max=project.eastingmax,
98#        northing_min=project.northingmin,
99#        northing_max= project.northingmax,
100        use_cache=True, 
101        verbose=True)
102
103#creates pts file for islands DEM
104dem2pts(islands_dir_name, use_cache=True, verbose=True)
105
[3851]106print'create Geospatial data objects from topographies'
[3827]107G1 = Geospatial_data(file_name = project.onshore_dir_name + '.pts')
108G2 = Geospatial_data(file_name = project.coast_dir_name + '.xya')
109G3 = Geospatial_data(file_name = project.islands_dir_name + '.pts')
110G_off = Geospatial_data(file_name = project.offshore_dir_name + '.xya')
111G_off1 = Geospatial_data(file_name = project.offshore_dir_name1 + '.xya')
112G_off2 = Geospatial_data(file_name = project.offshore_dir_name2 + '.xya')
113G_off3 = Geospatial_data(file_name = project.offshore_dir_name3 + '.xya')
114G_off4 = Geospatial_data(file_name = project.offshore_dir_name4 + '.xya')
115G_off5 = Geospatial_data(file_name = project.offshore_dir_name5 + '.xya')
116G_off6 = Geospatial_data(file_name = project.offshore_dir_name6 + '.xya')
117G_off7 = Geospatial_data(file_name = project.offshore_dir_name7 + '.xya')
118G_off8 = Geospatial_data(file_name = project.offshore_dir_name8 + '.xya')
119G_off9 = Geospatial_data(file_name = project.offshore_dir_name9 + '.xya')
120G_off10 = Geospatial_data(file_name = project.offshore_dir_name10 + '.xya')
121G_off11 = Geospatial_data(file_name = project.offshore_dir_name11 + '.xya')
122G_off12 = Geospatial_data(file_name = project.offshore_dir_name12 + '.xya')
123G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya')
124G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya')
125
[3851]126print'add all geospatial objects'
[3827]127G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \
128    + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \
129    + G_off13 + G_off14
130
[3851]131print'clip combined geospatial object by bounding polygon'
[3827]132G.clip(project.bounding_polygon)
133#FIXME: add a clip function to pts
134
[3851]135print'export combined DEM file'
136if access(project.topographies_time_dir,F_OK) == 0:
137    mkdir (project.topographies_time_dir)
138G.export_points_file(project.combined_time_dir_name + '.pts')
[3827]139
140#-------------------------------------------------------------------------
141# Convert URS to SWW file for boundary conditions
142#-------------------------------------------------------------------------
[3851]143print 'starting to create boundary conditions'
[3853]144boundaries_in_dir_name = project.boundaries_in_dir_name
[3827]145
[3851]146from anuga.shallow_water.data_manager import urs2sww
[3827]147
[3851]148print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
149print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
[3827]150
[3851]151if access(project.boundaries_time_dir,F_OK) == 0:
152    mkdir (project.boundaries_time_dir)
[3853]153urs2sww(boundaries_in_dir_name,basename_out= project.boundaries_time_dir_name,
[3851]154        minlat=project.south_boundary, maxlat=project.north_boundary,
[3827]155        minlon= project.west_boundary, maxlon=project.east_boundary,
[3851]156        mint=0, maxt= 35000,
[3827]157        verbose='true')
[3851]158       
[3827]159
160
161
162
163
164
[3851]165
166
Note: See TracBrowser for help on using the repository browser.