source: anuga_work/production/dampier_2008/build_boundary_27283.py @ 6026

Last change on this file since 6026 was 6026, checked in by kristy, 15 years ago

Had to run boundary script, have not set up the rest of the pythons scripts yet
File size: 7.8 KB
Line 
1"""
2Script for building boundary to run tsunami inundation scenario for busselton,
3WA, Australia.  The boundary is based on the National Hazard Map.
4
5Input: order_filename from project.py
6       event_number needs to be reflected in the project file
7Output: creates a sts file and csv files stored in project.boundaries_dir_event.
8The run_busselton.py is reliant on the output of this script.
9"""
10import os
11from os import sep
12
13
14from anuga.utilities.numerical_tools import ensure_numeric
15from Scientific.IO.NetCDF import NetCDFFile
16from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
17    compress,zeros,fabs,allclose,ones
18from anuga.utilities.polygon import inside_polygon,read_polygon
19from time import localtime, strftime, gmtime
20from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
21
22try:
23    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
24except:
25    print 'Cannot import pylab plotting will not work. Csv files are still created'
26
27#no project file
28#import project
29from os import getenv
30muxhome= getenv('MUXHOME')
31home = getenv('INUNDATIONHOME') + sep +'data'+sep
32state = 'western_australia'
33scenario_name = 'dampier'
34scenario = 'dampier_tsunami_scenario_2008'
35anuga_dir= home+state+sep+scenario+sep+'anuga'+sep
36boundaries_dir = anuga_dir+'boundaries'+ sep
37order_filename_dir = boundaries_dir + 'thinned_boundary_ordering.csv'
38boundaries_dir_event = boundaries_dir + '27283' + sep
39tide = 0
40#--------------------------------------------------------------------------
41# Create sts boundary from mux2files
42#--------------------------------------------------------------------------
43
44dir=os.path.join(muxhome,'mux')
45
46# Refer to event_027255.list in home+state+sep+event08 for event details
47# taken from David's event list for 27255
48urs_filenames = [
49    os.path.join(dir,'Java-0000-z.grd'),
50    os.path.join(dir,'Java-0001-z.grd'),
51    os.path.join(dir,'Java-0002-z.grd'),
52    os.path.join(dir,'Java-0003-z.grd'),
53    os.path.join(dir,'Java-0005-z.grd'), 
54    os.path.join(dir,'Java-0006-z.grd'), 
55    os.path.join(dir,'Java-0007-z.grd'), 
56    os.path.join(dir,'Java-0008-z.grd'), 
57    os.path.join(dir,'Java-0010-z.grd'), 
58    os.path.join(dir,'Java-0011-z.grd'), 
59    os.path.join(dir,'Java-0012-z.grd'), 
60    os.path.join(dir,'Java-0013-z.grd'), 
61    os.path.join(dir,'Java-0015-z.grd'), 
62    os.path.join(dir,'Java-0016-z.grd'), 
63    os.path.join(dir,'Java-0017-z.grd'), 
64    os.path.join(dir,'Java-0018-z.grd'), 
65    os.path.join(dir,'Java-0020-z.grd'), 
66    os.path.join(dir,'Java-0021-z.grd'), 
67    os.path.join(dir,'Java-0022-z.grd'), 
68    os.path.join(dir,'Java-0023-z.grd'), 
69    os.path.join(dir,'Java-0025-z.grd'), 
70    os.path.join(dir,'Java-0026-z.grd'), 
71    os.path.join(dir,'Java-0027-z.grd'), 
72    os.path.join(dir,'Java-0028-z.grd'), 
73    os.path.join(dir,'Java-0029-z.grd'), 
74    os.path.join(dir,'Java-0030-z.grd'), 
75    os.path.join(dir,'Java-0031-z.grd'), 
76    os.path.join(dir,'Java-0032-z.grd'), 
77    os.path.join(dir,'Java-0033-z.grd'), 
78    os.path.join(dir,'Java-0034-z.grd'), 
79    os.path.join(dir,'Java-0035-z.grd'), 
80    os.path.join(dir,'Java-0036-z.grd'), 
81    os.path.join(dir,'Java-0037-z.grd'), 
82    os.path.join(dir,'Java-0038-z.grd'), 
83    os.path.join(dir,'Java-0039-z.grd'), 
84    os.path.join(dir,'Java-0040-z.grd'), 
85    os.path.join(dir,'Java-0041-z.grd'), 
86    os.path.join(dir,'Java-0042-z.grd'), 
87    os.path.join(dir,'Java-0043-z.grd'), 
88    os.path.join(dir,'Java-0044-z.grd'), 
89    os.path.join(dir,'Java-0045-z.grd'), 
90    os.path.join(dir,'Java-0046-z.grd'), 
91    os.path.join(dir,'Java-0047-z.grd'),
92    os.path.join(dir,'Java-0048-z.grd')]   
93 
94print 'number of sources', len(urs_filenames)
95
96if len(urs_filenames) == 44: #change per event
97
98# AS per David Burbidge email on Friday 4th July the mag 9.3 event
99# has 1m worth of slip on each sub fault therefore mutliple each unit
100# source by the slip (10.4544) and sum the 44 time series together to
101# get the time series for this event at the points on your boundary.
102
103    weight_factor = 10.4544 #change per event
104    weights=weight_factor*ones(len(urs_filenames),Float)
105
106    #scenario_name=project.scenario_name 
107    order_filename=os.path.join(order_filename_dir) #project.
108    #boundaries_dir_event =project.boundaries_dir_event
109
110    print 'reading', order_filename
111    # Create ordered sts file
112    print 'creating sts file'
113
114    urs2sts(urs_filenames,
115            basename_out=os.path.join(boundaries_dir_event,scenario_name), 
116            ordering_filename=order_filename,
117            weights=weights,
118            mean_stage=tide,
119            verbose=True)
120else:
121    print 'number of sources do not match event.list'
122
123#------------------------------------------------------------------------------
124# Get gauges (timeseries of index points)
125#------------------------------------------------------------------------------
126def get_sts_gauge_data(filename,verbose=False):
127    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
128        compress,zeros,fabs,take
129    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
130    permutation = fid.variables['permutation'][:]
131    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
132    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
133    points=transpose(asarray([x.tolist(),y.tolist()]))
134    time=fid.variables['time'][:]+fid.starttime
135    elevation=fid.variables['elevation'][:]
136       
137    basename='sts_gauge'
138    quantity_names=['stage','xmomentum','ymomentum']
139    quantities = {}
140    for i, name in enumerate(quantity_names):
141        quantities[name] = fid.variables[name][:]
142
143#------------------------------------------------------------------------------
144# Get Maxium wave height throughout timeseries at each index point
145#------------------------------------------------------------------------------
146
147    maxname = 'max_sts_stage.csv'
148    fid_max = open(boundaries_dir_event+sep+maxname,'w') 
149    s = 'index, x, y, max_stage \n'
150    fid_max.write(s)   
151    for j in range(len(x)):
152        index= permutation[j]
153        stage = quantities['stage'][:,j]
154        xmomentum = quantities['xmomentum'][:,j]
155        ymomentum = quantities['ymomentum'][:,j]
156
157        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
158        fid_max.write(s)
159     
160#------------------------------------------------------------------------------
161# Get Minium wave height throughout timeseries at each index point
162#------------------------------------------------------------------------------
163
164    minname = 'min_sts_stage.csv'
165    fid_min = open(boundaries_dir_event+sep+minname,'w')
166    s = 'index, x, y, max_stage \n'
167    fid_min.write(s)   
168    for j in range(len(x)):
169        index= permutation[j]
170        stage = quantities['stage'][:,j]
171        xmomentum = quantities['xmomentum'][:,j]
172        ymomentum = quantities['ymomentum'][:,j]
173
174        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
175        fid_min.write(s)
176       
177
178
179        fid_sts = open(boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
180        s = 'time, stage, xmomentum, ymomentum \n'
181        fid_sts.write(s)
182
183#------------------------------------------------------------------------------
184# End of the get gauges
185#------------------------------------------------------------------------------
186        for k in range(len(time)-1):
187            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
188            fid_sts.write(s)
189
190        fid_sts.close()     
191
192    fid.close()
193
194
195    return quantities,elevation,time
196
197quantities,elevation,time=get_sts_gauge_data(os.path.join(boundaries_dir_event,scenario_name),verbose=False)
198
199print len(elevation), len(quantities['stage'][0,:])
200
Note: See TracBrowser for help on using the repository browser.