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

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

Current Hinwood scenario

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