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

Last change on this file since 5205 was 5205, checked in by herve, 16 years ago
File size: 5.2 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_clip.py
5The output sww file is stored in project_clip.output_time_dir
6
7The scenario is defined by a triangular mesh created from project_clip.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_clip   # 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_clip.output_build_time_dir,__file__, 
44               dirname(project_clip.__file__)+sep+ project_clip.__name__+'.py' )
45
46start_screen_catcher(project_clip.output_build_time_dir)
47
48print 'USER:    ', project_clip.user
49
50#-------------------------------------------------------------------------------
51# Preparation of topographic data# Convert ASC 2 DEM 2 PTS using source data and store result in source data
52# Do for coarse and fine data
53# Fine pts file to be clipped to area of interest
54#-------------------------------------------------------------------------------
55print"project_clip.combined_dir_name",project_clip.combined_dir_name
56
57# topography directory filenames
58total=project_clip.combined_dir_name
59
60print'create Geospatial data1 objects from topographies',total + '.txt'
61G1 = Geospatial_data(file_name = total + '.txt')
62
63print'add all geospatial objects'
64
65Gt= G1.clip_outside(Geospatial_data(project_clip.poly_west_reef)).clip_outside(Geospatial_data(project_clip.poly_west)).clip_outside(Geospatial_data(project_clip.poly_centerKolonia)).clip_outside(Geospatial_data(project_clip.poly_east))
66                         
67#print'clip combined geospatial object by bounding polygon'
68#G_clipped = G.clip(project_clip_urs.poly_all)
69#FIXME: add a clip function to pts
70#print'shape of clipped data', G_clipped.get_data_points().shape
71
72print'export combined DEM file'
73if access(project_clip.topographies_dir,F_OK) == 0:
74    mkdir (project_clip.topographies_dir)
75print'export',project_clip.combined_dir_name+ '.txt'
76Gt.export_points_file(project_clip.combined_dir_name1+ '.txt')
77
78
79
80'''
81print'project_clip.combined_dir_name + 1.xya',project_clip.combined_dir_name + '1.xya'
82G_all=Geospatial_data(file_name = project_clip.combined_dir_name + '1.xya')
83print'split'
84G_all_1, G_all_2 = G_all.split(.10)
85print'export 1'
86G_all_1.export_points_file(project_clip.combined_dir_name+'_small1' + '.xya')
87print'export 2'
88G_all_2.export_points_file(project_clip.combined_dir_name+'_other1' + '.xya')
89
90
91#-------------------------------------------------------------------------
92# Convert URS to SWW file for boundary conditions
93#-------------------------------------------------------------------------
94print 'starting to create boundary conditions'
95boundaries_in_dir_name = project_clip.boundaries_in_dir_name
96
97from anuga.shallow_water.data_manager import urs2sww
98
99print 'minlat=project_clip.north_boundary, maxlat=project_clip.south_boundary',project_clip.north_boundary, project_clip.south_boundary
100print 'minlon= project_clip.west_boundary, maxlon=project_clip.east_boundary',project_clip.west_boundary, project_clip.east_boundary
101
102#import sys; sys.exit()
103
104#if access(project_clip.boundaries_dir,F_OK) == 0:
105#    mkdir (project_clip.boundaries_dir)
106
107from caching import cache
108cache(urs2sww,
109      (boundaries_in_dir_name,
110       project_clip.boundaries_dir_name1),
111      {'verbose': True,
112       'minlat': project_clip.south_boundary,
113       'maxlat': project_clip.north_boundary,
114       'minlon': project_clip.west_boundary,
115       'maxlon': project_clip.east_boundary,
116#       'minlat': project_clip.south,
117#       'maxlat': project_clip.north,
118#       'minlon': project_clip.west,
119#       'maxlon': project_clip.east,
120       'mint': 0, 'maxt': 40000,
121#       'origin': domain.geo_reference.get_origin(),
122       'mean_stage': project_clip.tide,
123#       'zscale': 1,                 #Enhance tsunami
124       'fail_on_NaN': False},
125       verbose = True,
126       )
127#       dependencies = source_dir + project_clip.boundary_basename + '.sww')
128
129'''
130
131
132
133
134
135
136
137
Note: See TracBrowser for help on using the repository browser.