source: anuga_work/production/FSM/build_kolonia.py @ 5205

Last change on this file since 5205 was 5205, checked in by herve, 16 years ago
File size: 6.0 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.combined_dir_name",project.combined_dir_name
58
59# topography directory filenames
60onshore_dir_name = project.onshore_dir_name
61onshore_dir_name1 = project.onshore_dir_name1
62#island_in_dir_name = project_urs.island_in_dir_name
63Multibeam_dir_name = project.Multibeam_dir_name
64added_data_dir_name =project.added_data_dir_name
65
66print'create Geospatial data1 objects from topographies',onshore_dir_name + '.txt'
67G1 = Geospatial_data(file_name = onshore_dir_name + '.txt')
68
69print'create Geospatial data1 objects from topographies',onshore_dir_name1 + '.txt'
70G2 = Geospatial_data(file_name = onshore_dir_name1 + '.txt')
71
72print'create Geospatial data2 objects from multibeam', Multibeam_dir_name + '.txt'
73G_off1 = Geospatial_data(file_name = Multibeam_dir_name + '.txt')
74
75print'create Geospatial data objects from added data',added_data_dir_name + '.txt'
76G_off2 = Geospatial_data(file_name = added_data_dir_name + '.txt')
77
78print'add all geospatial objects'
79Gt1 = G1.clip(Geospatial_data(project.poly_west_reef))+G1.clip(Geospatial_data(project.poly_west))+G1.clip(Geospatial_data(project.poly_centerKolonia))+G1.clip(Geospatial_data(project.poly_east))
80Gt2= G2.clip(Geospatial_data(project.poly_topo))
81Gt3= G_off1.clip_outside(Geospatial_data(project.poly_west_reef)).clip_outside(Geospatial_data(project.poly_west)).clip_outside(Geospatial_data(project.poly_centerKolonia)).clip_outside(Geospatial_data(project.poly_east))
82                         
83G=Gt1+Gt2+Gt3+G_off2
84                         
85#print'clip combined geospatial object by bounding polygon'
86#G_clipped = G.clip(project_urs.poly_all)
87#FIXME: add a clip function to pts
88#print'shape of clipped data', G_clipped.get_data_points().shape
89
90print'export combined DEM file'
91if access(project.topographies_dir,F_OK) == 0:
92    mkdir (project.topographies_dir)
93print'export',project.combined_dir_name+ '.txt'
94G.export_points_file(project.combined_dir_name+ '.txt')
95
96
97
98'''
99print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
100G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
101print'split'
102G_all_1, G_all_2 = G_all.split(.10)
103print'export 1'
104G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
105print'export 2'
106G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
107
108
109#-------------------------------------------------------------------------
110# Convert URS to SWW file for boundary conditions
111#-------------------------------------------------------------------------
112print 'starting to create boundary conditions'
113boundaries_in_dir_name = project.boundaries_in_dir_name
114
115from anuga.shallow_water.data_manager import urs2sww
116
117print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
118print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
119
120#import sys; sys.exit()
121
122#if access(project.boundaries_dir,F_OK) == 0:
123#    mkdir (project.boundaries_dir)
124
125from caching import cache
126cache(urs2sww,
127      (boundaries_in_dir_name,
128       project.boundaries_dir_name1),
129      {'verbose': True,
130       'minlat': project.south_boundary,
131       'maxlat': project.north_boundary,
132       'minlon': project.west_boundary,
133       'maxlon': project.east_boundary,
134#       'minlat': project.south,
135#       'maxlat': project.north,
136#       'minlon': project.west,
137#       'maxlon': project.east,
138       'mint': 0, 'maxt': 40000,
139#       'origin': domain.geo_reference.get_origin(),
140       'mean_stage': project.tide,
141#       'zscale': 1,                 #Enhance tsunami
142       'fail_on_NaN': False},
143       verbose = True,
144       )
145#       dependencies = source_dir + project.boundary_basename + '.sww')
146
147'''
148
149
150
151
152
153
154
Note: See TracBrowser for help on using the repository browser.