source: anuga_work/production/australia_ph2/ceduna/build_boundary_27255.py @ 6291

Last change on this file since 6291 was 6291, checked in by myall, 14 years ago

new area for phase 2. scripts copied from standardised version

File size: 5.8 KB
Line 
1"""
2Script for building boundary to run a tsunami inundation scenario for busselton,
3WA, Australia.  The boundary is sourced from 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
12from os.path import join
13import project
14
15from anuga.utilities.numerical_tools import ensure_numeric
16from Scientific.IO.NetCDF import NetCDFFile
17from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
18    compress,zeros,fabs,allclose,ones
19from anuga.utilities.polygon import inside_polygon,read_polygon
20from time import localtime, strftime, gmtime
21from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
22
23try:
24    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
25except:
26    print 'Cannot import pylab plotting will not work. Csv files are still created'
27
28#--------------------------------------------------------------------------
29# Create sts boundary from mux2files
30#--------------------------------------------------------------------------
31# location of mux files
32# Note, this may change when all mux files are on the GA network
33dir=join(project.muxhome,'mux')
34
35# Refer to event_027255.list in home+state+sep+event08 for event details
36# taken from David's event list for 27255
37# Note, this is specific to the event and will be supplied by David.
38# In time, this will become more automated.
39urs_filenames = [
40    os.path.join(dir,'Java-0017-z.grd'), 
41    os.path.join(dir,'Java-0018-z.grd'),
42    os.path.join(dir,'Java-0019-z.grd'),
43    os.path.join(dir,'Java-0022-z.grd'), 
44    os.path.join(dir,'Java-0023-z.grd'), 
45    os.path.join(dir,'Java-0024-z.grd'), 
46    os.path.join(dir,'Java-0027-z.grd'), 
47    os.path.join(dir,'Java-0028-z.grd'), 
48    os.path.join(dir,'Java-0031-z.grd'), 
49    os.path.join(dir,'Java-0032-z.grd')] 
50 
51print 'number of sources', len(urs_filenames)
52
53if len(urs_filenames) == 10: #change per event
54
55    weight_factor = 30.84050 #change per event
56    weights=weight_factor*ones(len(urs_filenames),Float)
57
58    scenario_name = project.scenario_name
59    urs_order = join(project.urs_order)
60
61    print 'reading', urs_order
62
63    urs2sts(urs_filenames,
64            basename_out = join(project.event_sts),
65            ordering_filename = urs_order,
66            weights = weights,
67            mean_stage = project.tide,
68            verbose = True)
69else:
70    print 'number of sources do not match event.list'
71
72#------------------------------------------------------------------------------
73# Get gauges (timeseries of index points)
74#------------------------------------------------------------------------------
75def get_sts_gauge_data(filename,verbose=False):
76    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
77        compress,zeros,fabs,take
78    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
79    permutation = fid.variables['permutation'][:]
80    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
81    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
82    points=transpose(asarray([x.tolist(),y.tolist()]))
83    time=fid.variables['time'][:]+fid.starttime
84    elevation=fid.variables['elevation'][:]
85       
86    basename='sts_gauge'
87    quantity_names=['stage','xmomentum','ymomentum']
88    quantities = {}
89    for i, name in enumerate(quantity_names):
90        quantities[name] = fid.variables[name][:]
91
92    #------------------------------------------------------------------------------
93    # Get maxium wave height throughout timeseries at each index point
94    #------------------------------------------------------------------------------
95
96    maxname = 'max_sts_stage.csv'
97    fid_max = open(project.boundaries_dir_event+sep+maxname,'w')
98    s = 'index, x, y, max_stage \n'
99    fid_max.write(s)   
100    for j in range(len(x)):
101        index= permutation[j]
102        stage = quantities['stage'][:,j]
103        xmomentum = quantities['xmomentum'][:,j]
104        ymomentum = quantities['ymomentum'][:,j]
105
106        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
107        fid_max.write(s)
108     
109    #------------------------------------------------------------------------------
110    # Get minium wave height throughout timeseries at each index point
111    #------------------------------------------------------------------------------
112
113    minname = 'min_sts_stage.csv'
114    fid_min = open(project.boundaries_dir_event+sep+minname,'w')
115    s = 'index, x, y, max_stage \n'
116    fid_min.write(s)   
117    for j in range(len(x)):
118        index= permutation[j]
119        stage = quantities['stage'][:,j]
120        xmomentum = quantities['xmomentum'][:,j]
121        ymomentum = quantities['ymomentum'][:,j]
122
123        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
124        fid_min.write(s)
125       
126
127
128        fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
129        s = 'time, stage, xmomentum, ymomentum \n'
130        fid_sts.write(s)
131
132    #------------------------------------------------------------------------------
133    # End of the get gauges
134    #------------------------------------------------------------------------------
135        for k in range(len(time)-1):
136            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
137            fid_sts.write(s)
138
139        fid_sts.close()     
140
141    fid.close()
142
143
144    return quantities,elevation,time
145
146quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False)
147
148print len(elevation), len(quantities['stage'][0,:])
149
Note: See TracBrowser for help on using the repository browser.