1 | """ |
---|
2 | Replacement function that wraps the legacy function urs2sts() |
---|
3 | and passes in the data from the EventSelection <event>.list file. |
---|
4 | """ |
---|
5 | |
---|
6 | import os |
---|
7 | import os.path |
---|
8 | from time import localtime, strftime, gmtime |
---|
9 | |
---|
10 | import project |
---|
11 | |
---|
12 | from Scientific.IO.NetCDF import NetCDFFile |
---|
13 | import Numeric as num |
---|
14 | |
---|
15 | from anuga.shallow_water.data_manager import urs2sts |
---|
16 | |
---|
17 | import Numeric as num |
---|
18 | |
---|
19 | |
---|
20 | #------------------------------------------------------------------------------- |
---|
21 | # Get gauges (timeseries of index points) |
---|
22 | #------------------------------------------------------------------------------- |
---|
23 | def get_sts_gauge_data(filename, verbose=False): |
---|
24 | print 'get_sts_gauge_data: filename=%s' % filename |
---|
25 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
---|
26 | permutation = fid.variables['permutation'][:] |
---|
27 | x = fid.variables['x'][:] + fid.xllcorner #x-coordinates of vertices |
---|
28 | y = fid.variables['y'][:] + fid.yllcorner #y-coordinates of vertices |
---|
29 | points = num.transpose(num.asarray([x.tolist(), y.tolist()])) |
---|
30 | time = fid.variables['time'][:] + fid.starttime |
---|
31 | elevation = fid.variables['elevation'][:] |
---|
32 | |
---|
33 | basename = 'sts_gauge' |
---|
34 | quantity_names = ['stage', 'xmomentum', 'ymomentum'] |
---|
35 | quantities = {} |
---|
36 | for i, name in enumerate(quantity_names): |
---|
37 | quantities[name] = fid.variables[name][:] |
---|
38 | |
---|
39 | #--------------------------------------------------------------------------- |
---|
40 | # Get maximum wave height throughout timeseries at each index point |
---|
41 | #--------------------------------------------------------------------------- |
---|
42 | |
---|
43 | maxname = 'max_sts_stage.csv' |
---|
44 | print 'get_sts_gauge_data: maxname=%s' % maxname |
---|
45 | fid_max = open(os.path.join(project.event_folder, maxname), 'w') |
---|
46 | fid_max.write('index, x, y, max_stage \n') |
---|
47 | for j in range(len(x)): |
---|
48 | index = permutation[j] |
---|
49 | stage = quantities['stage'][:,j] |
---|
50 | xmomentum = quantities['xmomentum'][:,j] |
---|
51 | ymomentum = quantities['ymomentum'][:,j] |
---|
52 | |
---|
53 | fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage))) |
---|
54 | |
---|
55 | #--------------------------------------------------------------------------- |
---|
56 | # Get minimum wave height throughout timeseries at each index point |
---|
57 | #--------------------------------------------------------------------------- |
---|
58 | |
---|
59 | minname = 'min_sts_stage.csv' |
---|
60 | fid_min = open(os.path.join(project.event_folder, minname), 'w') |
---|
61 | fid_min.write('index, x, y, max_stage \n') |
---|
62 | for j in range(len(x)): |
---|
63 | index = permutation[j] |
---|
64 | stage = quantities['stage'][:,j] |
---|
65 | xmomentum = quantities['xmomentum'][:,j] |
---|
66 | ymomentum = quantities['ymomentum'][:,j] |
---|
67 | |
---|
68 | fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))) |
---|
69 | |
---|
70 | out_file = os.path.join(project.event_folder, |
---|
71 | basename+'_'+str(index)+'.csv') |
---|
72 | fid_sts = open(out_file, 'w') |
---|
73 | fid_sts.write('time, stage, xmomentum, ymomentum \n') |
---|
74 | |
---|
75 | #----------------------------------------------------------------------- |
---|
76 | # End of the get gauges |
---|
77 | #----------------------------------------------------------------------- |
---|
78 | for k in range(len(time)-1): |
---|
79 | fid_sts.write('%.6f, %.6f, %.6f, %.6f\n' |
---|
80 | % (time[k], stage[k], xmomentum[k], ymomentum[k])) |
---|
81 | |
---|
82 | fid_sts.close() |
---|
83 | |
---|
84 | fid.close() |
---|
85 | |
---|
86 | return quantities,elevation,time |
---|
87 | |
---|
88 | |
---|
89 | ## |
---|
90 | # @brief Build boundary STS files from one or more MUX files. |
---|
91 | # @param event_file Name of mux meta-file or single mux stem. |
---|
92 | # @param output_dir Directory to write STS data to. |
---|
93 | # @note 'event_file' is produced by EventSelection. |
---|
94 | def build_urs_boundary(event_file, output_dir): |
---|
95 | '''Build a boundary STS file from a set of MUX files.''' |
---|
96 | |
---|
97 | # if we are using an EventSelection multi-mux file |
---|
98 | if project.multi_mux: |
---|
99 | # get the mux+weight data from the meta-file (in <boundaries>) |
---|
100 | mux_event_file = os.path.join(project.event_folder, event_file) |
---|
101 | print 'using multi-mux file', mux_event_file |
---|
102 | try: |
---|
103 | fd = open(mux_event_file, 'r') |
---|
104 | mux_data = fd.readlines() |
---|
105 | fd.close() |
---|
106 | except IOError, e: |
---|
107 | msg = 'File %s cannot be read: %s' % (mux_event_file, str(e)) |
---|
108 | raise Exception, msg |
---|
109 | except: |
---|
110 | raise |
---|
111 | |
---|
112 | # first line of file is # filenames+weight in rest of file |
---|
113 | num_lines = int(mux_data[0].strip()) |
---|
114 | mux_data = mux_data[1:] |
---|
115 | print 'number of sources %d' % num_lines |
---|
116 | |
---|
117 | # quick sanity check on input mux meta-file |
---|
118 | if num_lines != len(mux_data): |
---|
119 | msg = ('Bad file %s: %d data lines, but line 1 count is %d' |
---|
120 | % (event_file, len(mux_data), num_lines)) |
---|
121 | raise Exception, msg |
---|
122 | |
---|
123 | # Create filename and weights lists. |
---|
124 | # Must chop GRD filename just after '*.grd'. |
---|
125 | mux_filenames = [] |
---|
126 | for line in mux_data: |
---|
127 | muxname = line.strip().split()[0] |
---|
128 | split_index = muxname.index('.grd') |
---|
129 | muxname = muxname[:split_index+len('.grd')] |
---|
130 | muxname = os.path.join(project.mux_data_folder, muxname) |
---|
131 | mux_filenames.append(muxname) |
---|
132 | |
---|
133 | mux_weights = [float(line.strip().split()[1]) for line in mux_data] |
---|
134 | |
---|
135 | # Call legacy function to create STS file. |
---|
136 | print 'creating sts file' |
---|
137 | urs2sts(mux_filenames, |
---|
138 | basename_out=output_dir, |
---|
139 | ordering_filename=project.urs_order, |
---|
140 | weights=mux_weights, |
---|
141 | central_meridian=project.central_meridian, |
---|
142 | zone=project.zone |
---|
143 | verbose=True) |
---|
144 | else: # a single mux stem file, assume 1.0 weight |
---|
145 | mux_file = os.path.join(project.event_folder, event_file) |
---|
146 | mux_filenames = [mux_file] |
---|
147 | print 'using single-mux file', mux_file |
---|
148 | |
---|
149 | weight_factor = 1.0 |
---|
150 | mux_weights = weight_factor*num.ones(len(mux_filenames), num.Float) |
---|
151 | |
---|
152 | order_filename = project.urs_order |
---|
153 | |
---|
154 | print 'reading', order_filename |
---|
155 | # Create ordered sts file |
---|
156 | urs2sts(mux_filenames, |
---|
157 | basename_out=output_dir, |
---|
158 | ordering_filename=order_filename, |
---|
159 | weights=mux_weights, |
---|
160 | central_meridian=project.central_meridian, |
---|
161 | zone=project.zone |
---|
162 | verbose=True) |
---|
163 | |
---|
164 | # report on progress so far |
---|
165 | sts_file = os.path.join(project.event_folder, project.scenario_name) |
---|
166 | quantities, elevation, time = get_sts_gauge_data(sts_file, verbose=False) |
---|
167 | print len(elevation), len(quantities['stage'][0,:]) |
---|
168 | |
---|
169 | |
---|
170 | |
---|