source: anuga_work/production/Broome_2008/build_boundary_27255.py @ 5810

Last change on this file since 5810 was 5810, checked in by rwilson, 15 years ago

Initial setup of the new Broome scenario.

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