source: trunk/anuga_validation/automated_validation_tests/patong_beach_validation/build_urs_boundary.py @ 8427

Last change on this file since 8427 was 8427, checked in by davies, 12 years ago

Adding the trapezoidal channel validation test, and editing the ANUGA manual

File size: 6.5 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 numpy as num
14
15from anuga.shallow_water.data_manager import urs2sts
16
17
18#-------------------------------------------------------------------------------
19# Get gauges (timeseries of index points)
20#-------------------------------------------------------------------------------
21def get_sts_gauge_data(filename, verbose=False):
22    print 'get_sts_gauge_data: filename=%s' % filename
23    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
24    permutation = fid.variables['permutation'][:]
25    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
26    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
27    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
28    time = fid.variables['time'][:] + fid.starttime
29    elevation = fid.variables['elevation'][:]
30       
31    basename = 'sts_gauge'
32    quantity_names = ['stage', 'xmomentum', 'ymomentum']
33    quantities = {}
34    for i, name in enumerate(quantity_names):
35        quantities[name] = fid.variables[name][:]
36
37    #---------------------------------------------------------------------------
38    # Get maximum wave height throughout timeseries at each index point
39    #---------------------------------------------------------------------------
40     
41    maxname = 'max_sts_stage.csv'
42    print 'get_sts_gauge_data: maxname=%s' % maxname
43    fid_max = open(os.path.join(project.event_folder, 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_folder, 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_folder,
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 event_file Name of mux meta-file or single mux stem.
90# @param output_dir Directory to write STS data to.
91# @note 'event_file' is produced by EventSelection.
92def build_urs_boundary(event_file, output_dir):
93    '''Build a boundary STS file from a set of MUX files.'''
94
95    # if we are using an EventSelection multi-mux file
96    if project.multi_mux:
97        # get the mux+weight data from the meta-file (in <boundaries>)
98        mux_event_file = os.path.join(project.event_folder, event_file)
99        print 'using multi-mux file', mux_event_file
100        try:
101            fd = open(mux_event_file, 'r')
102            mux_data = fd.readlines()
103            fd.close()
104        except IOError, e:
105            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
106            raise Exception, msg
107        except:
108            raise
109
110        # first line of file is # filenames+weight in rest of file
111        num_lines = int(mux_data[0].strip())
112        mux_data = mux_data[1:]
113        print 'number of sources %d' % num_lines
114
115        # quick sanity check on input mux meta-file
116        if num_lines != len(mux_data):
117            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
118                   % (event_file, len(mux_data), num_lines))
119            raise Exception, msg
120
121        # Create filename and weights lists.
122        # Must chop GRD filename just after '*.grd'.
123        mux_filenames = []
124        for line in mux_data:
125            muxname = line.strip().split()[0]
126            split_index = muxname.index('.grd')
127            muxname = muxname[:split_index+len('.grd')]
128            muxname = os.path.join(project.mux_data_folder, muxname)
129            mux_filenames.append(muxname)
130       
131        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
132
133        # Call legacy function to create STS file.
134        print 'creating sts file'
135        urs2sts(mux_filenames,
136                basename_out=output_dir,
137                ordering_filename=project.urs_order,
138                weights=mux_weights,
139                zone=project.zone,
140                mean_stage=project.tide,
141                verbose=True)
142    else:                           # a single mux stem file, assume 1.0 weight
143        mux_file = os.path.join(project.event_folder, event_file)
144        mux_filenames = [mux_file]
145        print 'using single-mux file', mux_file
146
147        weight_factor = 1.0
148        mux_weights = weight_factor*num.ones(len(mux_filenames), num.float)
149           
150        order_filename = project.urs_order
151
152        print 'reading', order_filename
153        # Create ordered sts file
154        urs2sts(mux_filenames,
155                basename_out=output_dir,
156                ordering_filename=order_filename,
157                weights=mux_weights,
158                mean_stage=project.tide,
159                verbose=True)
160
161    # report on progress so far
162    sts_file = os.path.join(project.event_folder, project.scenario_name)
163    quantities, elevation, time = get_sts_gauge_data(sts_file, verbose=False)
164    print len(elevation), len(quantities['stage'][0,:])
165
166   
167
Note: See TracBrowser for help on using the repository browser.