| 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 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
|---|
| 25 | permutation = fid.variables['permutation'][:] |
|---|
| 26 | x = fid.variables['x'][:] + fid.xllcorner #x-coordinates of vertices |
|---|
| 27 | y = fid.variables['y'][:] + fid.yllcorner #y-coordinates of vertices |
|---|
| 28 | points = num.transpose(num.asarray([x.tolist(), y.tolist()])) |
|---|
| 29 | time = fid.variables['time'][:] + fid.starttime |
|---|
| 30 | elevation = fid.variables['elevation'][:] |
|---|
| 31 | |
|---|
| 32 | basename = 'sts_gauge' |
|---|
| 33 | quantity_names = ['stage', 'xmomentum', 'ymomentum'] |
|---|
| 34 | quantities = {} |
|---|
| 35 | for i, name in enumerate(quantity_names): |
|---|
| 36 | quantities[name] = fid.variables[name][:] |
|---|
| 37 | |
|---|
| 38 | #--------------------------------------------------------------------------- |
|---|
| 39 | # Get maximum wave height throughout timeseries at each index point |
|---|
| 40 | #--------------------------------------------------------------------------- |
|---|
| 41 | |
|---|
| 42 | maxname = 'max_sts_stage.csv' |
|---|
| 43 | fid_max = open(os.path.join(project.event_sts, maxname), 'w') |
|---|
| 44 | fid_max.write('index, x, y, max_stage \n') |
|---|
| 45 | for j in range(len(x)): |
|---|
| 46 | index = permutation[j] |
|---|
| 47 | stage = quantities['stage'][:,j] |
|---|
| 48 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 49 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 50 | |
|---|
| 51 | fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage))) |
|---|
| 52 | |
|---|
| 53 | #--------------------------------------------------------------------------- |
|---|
| 54 | # Get minimum wave height throughout timeseries at each index point |
|---|
| 55 | #--------------------------------------------------------------------------- |
|---|
| 56 | |
|---|
| 57 | minname = 'min_sts_stage.csv' |
|---|
| 58 | fid_min = open(os.path.join(project.event_sts, minname), 'w') |
|---|
| 59 | fid_min.write('index, x, y, max_stage \n') |
|---|
| 60 | for j in range(len(x)): |
|---|
| 61 | index = permutation[j] |
|---|
| 62 | stage = quantities['stage'][:,j] |
|---|
| 63 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 64 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 65 | |
|---|
| 66 | fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))) |
|---|
| 67 | |
|---|
| 68 | out_file = os.path.join(project.event_sts, |
|---|
| 69 | basename+'_'+str(index)+'.csv') |
|---|
| 70 | fid_sts = open(out_file, 'w') |
|---|
| 71 | fid_sts.write('time, stage, xmomentum, ymomentum \n') |
|---|
| 72 | |
|---|
| 73 | #----------------------------------------------------------------------- |
|---|
| 74 | # End of the get gauges |
|---|
| 75 | #----------------------------------------------------------------------- |
|---|
| 76 | for k in range(len(time)-1): |
|---|
| 77 | fid_sts.write('%.6f, %.6f, %.6f, %.6f\n' |
|---|
| 78 | % (time[k], stage[k], xmomentum[k], ymomentum[k])) |
|---|
| 79 | |
|---|
| 80 | fid_sts.close() |
|---|
| 81 | |
|---|
| 82 | fid.close() |
|---|
| 83 | |
|---|
| 84 | return quantities,elevation,time |
|---|
| 85 | |
|---|
| 86 | |
|---|
| 87 | ## |
|---|
| 88 | # @brief Build boundary STS files from one or more MUX files. |
|---|
| 89 | # @param mux_dir Directory containing the MUX files. |
|---|
| 90 | # @param event_file Path to meta-file containing mux files+weight data. |
|---|
| 91 | # @param output_dir Directory to write STS data to. |
|---|
| 92 | # @note 'event_file' is produced by EventSelection. |
|---|
| 93 | def build_urs_boundary(mux_dir, event_file, output_dir): |
|---|
| 94 | '''Build a boundary STS file from a set of MUX files.''' |
|---|
| 95 | |
|---|
| 96 | print 'build_urs_boundary: mux_dir=%s' % mux_dir |
|---|
| 97 | print 'build_urs_boundary: event_file=%s' % event_file |
|---|
| 98 | print 'build_urs_boundary: output_dir=%s' % output_dir |
|---|
| 99 | |
|---|
| 100 | # if we are using an EventSelection multi-mux file |
|---|
| 101 | if project.multi_mux: |
|---|
| 102 | # get the mux+weight data from the event file |
|---|
| 103 | mux_event_file = os.path.join(mux_dir, 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(mux_dir, 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 | mux_filenames = [os.path.join(mux_dir, event_file)] |
|---|
| 147 | |
|---|
| 148 | weight_factor = 1.0 |
|---|
| 149 | mux_weights = weight_factor*num.ones(len(mux_filenames), num.Float) |
|---|
| 150 | |
|---|
| 151 | order_filename = project.order_filename_dir |
|---|
| 152 | |
|---|
| 153 | print 'reading', order_filename |
|---|
| 154 | # Create ordered sts file |
|---|
| 155 | print 'creating sts file' |
|---|
| 156 | urs2sts(mux_filenames, |
|---|
| 157 | basename_out=output_dir, |
|---|
| 158 | ordering_filename=order_filename, |
|---|
| 159 | weights=mux_weights, |
|---|
| 160 | mean_stage=project.tide, |
|---|
| 161 | verbose=True) |
|---|
| 162 | |
|---|
| 163 | # report on progress so far |
|---|
| 164 | quantities, elevation, time = get_sts_gauge_data(project.event_folder, |
|---|
| 165 | verbose=False) |
|---|
| 166 | print len(elevation), len(quantities['stage'][0,:]) |
|---|
| 167 | |
|---|
| 168 | |
|---|
| 169 | |
|---|
