source: anuga_work/production/busselton/build_boundary_27283.py @ 5786

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

Updated all scripts to reflect Perth format. Also added 250m scripts

File size: 7.4 KB
RevLine 
[5786]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
12import project
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#--------------------------------------------------------------------------
28# Create sts boundary from mux2files
29#--------------------------------------------------------------------------
30
31dir=os.path.join(project.muxhome,'mux')
32
33# Refer to event_027255.list in home+state+sep+event08 for event details
34# taken from David's event list for 27255
35urs_filenames = [
36    os.path.join(dir,'Java-0000-z.grd'),
37    os.path.join(dir,'Java-0001-z.grd'),
38    os.path.join(dir,'Java-0002-z.grd'),
39    os.path.join(dir,'Java-0003-z.grd'),
40    os.path.join(dir,'Java-0005-z.grd'), 
41    os.path.join(dir,'Java-0006-z.grd'), 
42    os.path.join(dir,'Java-0007-z.grd'), 
43    os.path.join(dir,'Java-0008-z.grd'), 
44    os.path.join(dir,'Java-0010-z.grd'), 
45    os.path.join(dir,'Java-0011-z.grd'), 
46    os.path.join(dir,'Java-0012-z.grd'), 
47    os.path.join(dir,'Java-0013-z.grd'), 
48    os.path.join(dir,'Java-0015-z.grd'), 
49    os.path.join(dir,'Java-0016-z.grd'), 
50    os.path.join(dir,'Java-0017-z.grd'), 
51    os.path.join(dir,'Java-0018-z.grd'), 
52    os.path.join(dir,'Java-0020-z.grd'), 
53    os.path.join(dir,'Java-0021-z.grd'), 
54    os.path.join(dir,'Java-0022-z.grd'), 
55    os.path.join(dir,'Java-0023-z.grd'), 
56    os.path.join(dir,'Java-0025-z.grd'), 
57    os.path.join(dir,'Java-0026-z.grd'), 
58    os.path.join(dir,'Java-0027-z.grd'), 
59    os.path.join(dir,'Java-0028-z.grd'), 
60    os.path.join(dir,'Java-0029-z.grd'), 
61    os.path.join(dir,'Java-0030-z.grd'), 
62    os.path.join(dir,'Java-0031-z.grd'), 
63    os.path.join(dir,'Java-0032-z.grd'), 
64    os.path.join(dir,'Java-0033-z.grd'), 
65    os.path.join(dir,'Java-0034-z.grd'), 
66    os.path.join(dir,'Java-0035-z.grd'), 
67    os.path.join(dir,'Java-0036-z.grd'), 
68    os.path.join(dir,'Java-0037-z.grd'), 
69    os.path.join(dir,'Java-0038-z.grd'), 
70    os.path.join(dir,'Java-0039-z.grd'), 
71    os.path.join(dir,'Java-0040-z.grd'), 
72    os.path.join(dir,'Java-0041-z.grd'), 
73    os.path.join(dir,'Java-0042-z.grd'), 
74    os.path.join(dir,'Java-0043-z.grd'), 
75    os.path.join(dir,'Java-0044-z.grd'), 
76    os.path.join(dir,'Java-0045-z.grd'), 
77    os.path.join(dir,'Java-0046-z.grd'), 
78    os.path.join(dir,'Java-0047-z.grd'),
79    os.path.join(dir,'Java-0048-z.grd')]   
80 
81print 'number of sources', len(urs_filenames)
82
83if len(urs_filnames) = 44 #change per event
84
85# AS per David Burbidge email on Friday 4th July the mag 9.3 event
86# has 1m worth of slip on each sub fault therefore mutliple each unit
87# source by the slip (10.4544) and sum the 44 time series together to
88# get the time series for this event at the points on your boundary.
89
90    weight_factor = 10.4544 #change per event
91    weights=weight_factor*ones(len(urs_filenames),Float)
92
93    scenario_name=project.scenario_name
94    order_filename=os.path.join(project.order_filename_dir)
95
96    print 'reading', order_filename
97    # Create ordered sts file
98    print 'creating sts file'
99
100    urs2sts(urs_filenames,
101            basename_out=os.path.join(project.boundaries_dir_event,scenario_name),
102            ordering_filename=order_filename,
103            weights=weights,
104            mean_stage=project.tide,
105            verbose=True)
106else print 'number of sources do not match event.list'
107
108#------------------------------------------------------------------------------
109# Get gauges (timeseries of index points)
110#------------------------------------------------------------------------------
111def get_sts_gauge_data(filename,verbose=False):
112    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
113        compress,zeros,fabs,take
114    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
115    permutation = fid.variables['permutation'][:]
116    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
117    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
118    points=transpose(asarray([x.tolist(),y.tolist()]))
119    time=fid.variables['time'][:]+fid.starttime
120    elevation=fid.variables['elevation'][:]
121       
122    basename='sts_gauge'
123    quantity_names=['stage','xmomentum','ymomentum']
124    quantities = {}
125    for i, name in enumerate(quantity_names):
126        quantities[name] = fid.variables[name][:]
127
128#------------------------------------------------------------------------------
129# Get Maxium wave height throughout timeseries at each index point
130#------------------------------------------------------------------------------
131
132    maxname = 'max_sts_stage.csv'
133    fid_max = open(project.boundaries_dir_event+sep+maxname,'w')
134    s = 'index, x, y, max_stage \n'
135    fid_max.write(s)   
136    for j in range(len(x)):
137        index= permutation[j]
138        stage = quantities['stage'][:,j]
139        xmomentum = quantities['xmomentum'][:,j]
140        ymomentum = quantities['ymomentum'][:,j]
141
142        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
143        fid_max.write(s)
144     
145#------------------------------------------------------------------------------
146# Get Minium wave height throughout timeseries at each index point
147#------------------------------------------------------------------------------
148
149    minname = 'min_sts_stage.csv'
150    fid_min = open(project.boundaries_dir_event+sep+maxname,'w')
151    s = 'index, x, y, max_stage \n'
152    fid_min.write(s)   
153    for j in range(len(x)):
154        index= permutation[j]
155        stage = quantities['stage'][:,j]
156        xmomentum = quantities['xmomentum'][:,j]
157        ymomentum = quantities['ymomentum'][:,j]
158
159        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
160        fid_min.write(s)
161       
162
163
164        fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
165        s = 'time, stage, xmomentum, ymomentum \n'
166        fid_sts.write(s)
167
168#------------------------------------------------------------------------------
169# End of the get gauges
170#------------------------------------------------------------------------------
171        for k in range(len(time)-1):
172            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
173            fid_sts.write(s)
174
175        fid_sts.close()     
176
177    fid.close()
178
179
180    return quantities,elevation,time
181
182quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False)
183
184print len(elevation), len(quantities['stage'][0,:])
185
Note: See TracBrowser for help on using the repository browser.