source: anuga_work/production/carnarvon/build_boundary_27283.py @ 5790

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

Updated scripts to reflect perth. Also added the 250m scripts and used more than one polygon for initial condition
File size: 7.4 KB
Line 
1"""
2Script for building boundary to run tsunami inundation scenario for carnarvon,
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_carnarvon.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.