source: anuga_work/development/Hinwood_2008/prepare_time_boundary.py @ 5455

Last change on this file since 5455 was 5455, checked in by duncan, 15 years ago

Current Hinwood scenario - cropping the flume length

File size: 8.4 KB
RevLine 
[5076]1"""
2Script for running a breaking wave simulation of Jon Hinwoods wave tank.
3Note: this is based on the frinction_ua_flume_2006 structure.
4
5
6Duncan Gray, GA - 2007
7
8"""
9
10
11#----------------------------------------------------------------------------
12# Import necessary modules
13#----------------------------------------------------------------------------
14
15from Scientific.IO.NetCDF import NetCDFFile
16from Numeric import array, zeros, Float
[5350]17 
18from anuga.utilities.numerical_tools import ensure_numeric
[5076]19
[5350]20from interp import interp
[5076]21
[5392]22from os.path import join
23
24# from os import getenv
25# from os.path import join
26# home = getenv('INUNDATIONHOME')
27# Hinwood_dir = join(home,'data','flumes','Hinwood_2008')
28# raw_data_dir = join(Hinwood_dir, 'raw_data')
29   
30def csv2tms(filename, offshore_bed_elevation):
[5076]31    """Convert benchmark 2 time series to NetCDF tms file.
[5392]32    the filename is the name of the output tms file, eg 'hi.tsm'.
33    There must be an equivalent .csv file, eg 'hi.csv'.
[5370]34   
[5076]35    """
36
37
38
39    print 'Creating', filename
40
41    # Read the ascii (.csv) version of this file
42    fid = open(filename[:-4] + '.csv')
43
44    # Read remaining lines
45    lines = fid.readlines()
46    fid.close()
47
48
49    N = len(lines)
50    T = zeros(N, Float)  #Time
51    Q = zeros(N, Float)  #Stage
52    X = zeros(N, Float)  #XMomentum
53    Y = zeros(N, Float)  #YMomentum
54
55    for i, line in enumerate(lines):
56        fields = line.split(',')
57
58        T[i] = float(fields[0])
[5392]59        Q[i] = float(fields[1])
60        depth = Q[i] - offshore_bed_elevation
[5076]61        X[i] = float(fields[2]) * depth
[5350]62        try:
63            Y[i] = float(fields[3]) * depth
64        except:
65            pass
[5076]66
67
[5370]68   
[5076]69    # Create tms NetCDF file
70    fid = NetCDFFile(filename, 'w')
71    fid.institution = 'Geoscience Australia'
72    fid.description = 'Input wave for Benchmark 2'
73    fid.starttime = 0.0
74    fid.createDimension('number_of_timesteps', len(T))
75    fid.createVariable('time', Float, ('number_of_timesteps',))
76    fid.variables['time'][:] = T
77
78    fid.createVariable('stage', Float, ('number_of_timesteps',))
79    fid.variables['stage'][:] = Q[:]
80
81    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
82    fid.variables['xmomentum'][:] =  X[:]
83
84    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
85    fid.variables['ymomentum'][:] =  Y[:]
86
87    fid.close()
88
[5350]89
[5395]90def combine_velocity_depth(velocity_file, depth_file, out_file, metadata_dic):
[5350]91    """
92   
[5370]93    Convert the raw velocity and depth values, which have values at
[5350]94    different times to a csv file, with values at the same time, with
95    SI units.
96
[5370]97    Set the depth values to be at the same times as the velocity values.
98
[5350]99    The format for the velocity file is;
[5395]100    [time, sec], [x-velocity +ve is towards the wave generator, cm/sec],
[5350]101    [y-velocity], [z-velocity]
102
103    The format for the pressure file is
104    [time, sec], [mm above SWL for sensor A], many other sensors...
[5395]105
106    The format of the output file is;
107    [time, sec], [m above SWL for sensor A],
108    [x-velocity +ve is towards the shore, m/sec],
109    [y-velocity +ve is towards the left of the tank, looking from the
110    generator to the shore, m/sec]
111   
112   
[5350]113    """
[5370]114    missing = 1e+20
115   
[5350]116    # Read velocity file
117    vfid = open(velocity_file)
118    lines = vfid.readlines()
119    vfid.close()
120
121   
122    n_velocity = len(lines)
123    vtimes = zeros(n_velocity, Float)  #Time
[5370]124    x_velocities = zeros(n_velocity, Float)  #
125    y_velocities = zeros(n_velocity, Float)  #
[5350]126
127    for i, line in enumerate(lines):
[5370]128        fields = line.split() #(',')
[5350]129
130        vtimes[i] = float(fields[0])
[5370]131        x_velocities[i] = float(fields[1])
132        y_velocities[i] = float(fields[2])
[5350]133
134    # Read the depth file
135    dfid = open(depth_file)
136    lines = dfid.readlines()
137    dfid.close()
138
[5370]139   
[5350]140    n_depth = len(lines)
[5370]141    n_sensors = len(lines[0].split())
[5350]142    dtimes = zeros(n_depth, Float)  #Time
143    depths = zeros(n_depth, Float)  #
[5370]144    sensors = zeros((n_depth,n_sensors), Float)
[5350]145   
146    for i, line in enumerate(lines):
[5370]147        fields = line.split() #(',')
148        fields = [float(j) for j in fields]
149        dtimes[i] = fields[0]
150        depths[i] = fields[1]
151        sensors[i] = fields
152    #print "dtimes", dtimes
153    #print "depths", depths
154    #print "vtimes", vtimes
155    depths_at_vtimes = interp( depths, dtimes,
156                               vtimes, missing=missing)
157    depths_at_vtimes = ensure_numeric(depths_at_vtimes)
[5350]158
[5370]159    #print "len(dtimes)", len(vtimes)
160    #print "len(depths_at_vtimes)", len(depths_at_vtimes)
161#     for i in range(len(depths_at_vtimes)):
162#         print "i", i
163#         print "vtimes[i]", vtimes[i]
164#         print "depths_at_vtimes[i]", depths_at_vtimes[i]
165#     print "depths_at_vtimes",  depths_at_vtimes
[5350]166    depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m
[5370]167    missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes
168    x_velocities = ensure_numeric(x_velocities)
[5395]169    # Swap axis around convert cm/sec to m/sec
170    x_velocities = x_velocities * -0.01 
[5370]171    y_velocities = ensure_numeric(y_velocities)
[5395]172    # Swap axis around convert cm/sec to m/sec
173    y_velocities = y_velocities * -0.01
[5350]174   
175    fid = open(out_file,'w')
176
177    assert len(depths_at_vtimes) == len(vtimes)
[5370]178    start_time = None
179    #start_time = 0.0
[5350]180    #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes,
181    #                                           velocities):
182    for i in xrange(len(vtimes)):
[5370]183        if not depths_at_vtimes[i] == missing:
184            # Make the times start at zero.
185            if start_time is None:
[5395]186                start_time = vtimes[i]
187            final_time = vtimes[i]-start_time
188            fid.write(str(final_time) \
[5370]189                      + ',' + str(depths_at_vtimes[i]) \
190                      + ',' + str(x_velocities[i]) \
191                      + ',' + str(y_velocities[i])+'\n')
[5350]192
193    fid.close()
194   
[5370]195    # Since there is a new time reference save the depth info using this
196    # new reference.
[5395]197    fid = open(depth_file[:-4] + '_exp_depth.csv','w')
[5370]198    sensors[:,0] -= start_time
199    #print "depth_file", depth_file
200    #print "start_time", start_time
[5395]201
[5455]202   
203    bed_elevation_list = metadata_dic['gauge_bed_elevation']
204    # +2 due to Time column and gauge A
205    max_j = len(bed_elevation_list)+2
[5395]206    # Write a header
207    fid.write('Time')
208    for gauge_x in metadata_dic['gauge_x']:
209        fid.write(',' + str(gauge_x)) 
210    fid.write('\n')
211    for i in xrange(len(dtimes)):       
212        fid.write(str(sensors[i,0])) # Time
213        # Don't write sensor A.
214        # It is the boundary condition.
[5455]215        for j, bed_elevation in map(None,
216                                    xrange(2,max_j),
217                                    bed_elevation_list):
[5395]218            # depth, m
219            gauge = sensors[i,j]/1000 - bed_elevation # Convert from mm to m
220            fid.write(',' + str(gauge)) 
[5370]221        fid.write('\n')
222    fid.close()
[5455]223
224   
225    # Since there is a new time reference save the stage info using this
226    # new reference.
227    fid = open(depth_file[:-4] + '_exp_stage.csv','w')
228    sensors[:,0] -= start_time
229   
230    # Write a header
231    fid.write('Time')
232    for gauge_x in metadata_dic['gauge_x']:
233        fid.write(',' + str(gauge_x)) 
234    fid.write('\n')
235    for i in xrange(len(dtimes)):       
236        fid.write(str(sensors[i,0])) # Time
237        # Don't write sensor A.
238        # It is the boundary condition.
239        for j in xrange(2,max_j):
240            # stage, m
241            gauge = sensors[i,j]/1000 # Convert from mm to m
242            fid.write(',' + str(gauge)) 
243        fid.write('\n')
244    fid.close()
245   
[5395]246    return final_time
[5370]247
[5392]248def prepare_time_boundary(metadata_dic, raw_data_dir, output_dir):
249    """
250    """
251    scenario_id = metadata_dic['scenario_id']
252    velocity_file = join(raw_data_dir,scenario_id+'velfilt.txt')
253    depth_file = join(raw_data_dir,scenario_id+'pressfilt.txt')
254    out_file = join(output_dir, scenario_id+'_boundary.csv')
[5395]255    #out_depth_file = join(output_dir, scenario_id+'_exp_depth.csv')
[5350]256   
[5395]257    final_time = combine_velocity_depth(velocity_file, depth_file, out_file,
258                                        metadata_dic) 
[5392]259    if metadata_dic['xleft'][1] >= 0.0:
260        # This should be a -ve value, since the still water level is the
261        # z origin.
[5395]262        print "Warning: The z origin seems incorrect."     
263    tsm_file = out_file[:-4] + '.tsm' 
[5392]264    csv2tms(tsm_file, metadata_dic['xleft'][1])
[5395]265    return final_time
[5370]266   
[5350]267#-------------------------------------------------------------------
[5076]268if __name__ == "__main__":
[5392]269    pass
Note: See TracBrowser for help on using the repository browser.