source: anuga_work/production/busselton/standardised_version/build_urs_boundary.py @ 6329

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

Further changes from testing.

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