source: anuga_work/production/perth/build_boundary_27255.py @ 5804

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