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

"""

Script for running a breaking wave simulation of Jon Hinwoods wave tank.
Note: this is based on the frinction_ua_flume_2006 structure.


Duncan Gray, GA - 2007 



"""


#----------------------------------------------------------------------------
# Import necessary modules
#----------------------------------------------------------------------------

from Scientific.IO.NetCDF import NetCDFFile
from Numeric import array, zeros, Float
  
from anuga.utilities.numerical_tools import ensure_numeric

from interp import interp

def prepare_time_boundary(filename):
    """Convert benchmark 2 time series to NetCDF tms file.
   
    """



    print 'Creating', filename

    # Read the ascii (.csv) version of this file
    fid = open(filename[:-4] + '.csv')

    # Read remaining lines
    lines = fid.readlines()
    fid.close()


    N = len(lines)
    T = zeros(N, Float)  #Time
    Q = zeros(N, Float)  #Stage
    X = zeros(N, Float)  #XMomentum
    Y = zeros(N, Float)  #YMomentum

    for i, line in enumerate(lines):
        fields = line.split(',')

        T[i] = float(fields[0])
        Q[i] = depth = float(fields[1])
        X[i] = float(fields[2]) * depth
        try:
            Y[i] = float(fields[3]) * depth
        except:
            pass


    
    # Create tms NetCDF file
    fid = NetCDFFile(filename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for Benchmark 2'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', Float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', Float, ('number_of_timesteps',))
    fid.variables['stage'][:] = Q[:]

    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] =  X[:]

    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
    fid.variables['ymomentum'][:] =  Y[:]

    fid.close()


def combine_velocity_depth(velocity_file, depth_file, out_file):
    """
    
    Convert the raw velocity and depth values, which have values at
    different times to a csv file, with values at the same time, with
    SI units.

    Set the depth values to be at the same times as the velocity values.

    The format for the velocity file is;
    [time, sec], [x-velocity +ve is towards the wave generator, m/sec],
    [y-velocity], [z-velocity]

    The format for the pressure file is
    [time, sec], [mm above SWL for sensor A], many other sensors...
    """
    missing = 1e+20
    
    # Read velocity file
    vfid = open(velocity_file)
    lines = vfid.readlines()
    vfid.close()

    
    n_velocity = len(lines)
    vtimes = zeros(n_velocity, Float)  #Time
    x_velocities = zeros(n_velocity, Float)  #
    y_velocities = zeros(n_velocity, Float)  #

    for i, line in enumerate(lines):
        fields = line.split() #(',')

        vtimes[i] = float(fields[0])
        x_velocities[i] = float(fields[1])
        y_velocities[i] = float(fields[2])

    # Read the depth file
    dfid = open(depth_file)
    lines = dfid.readlines()
    dfid.close()

    
    n_depth = len(lines)
    n_sensors = len(lines[0].split())
    dtimes = zeros(n_depth, Float)  #Time
    depths = zeros(n_depth, Float)  #
    sensors = zeros((n_depth,n_sensors), Float)
    
    for i, line in enumerate(lines):
        fields = line.split() #(',')
        fields = [float(j) for j in fields]
        dtimes[i] = fields[0]
        depths[i] = fields[1]
        sensors[i] = fields
    #print "dtimes", dtimes
    #print "depths", depths
    #print "vtimes", vtimes
    depths_at_vtimes = interp( depths, dtimes,
                               vtimes, missing=missing)
    depths_at_vtimes = ensure_numeric(depths_at_vtimes)

    #print "len(dtimes)", len(vtimes)
    #print "len(depths_at_vtimes)", len(depths_at_vtimes)
#     for i in range(len(depths_at_vtimes)):
#         print "i", i
#         print "vtimes[i]", vtimes[i]
#         print "depths_at_vtimes[i]", depths_at_vtimes[i]
#     print "depths_at_vtimes",  depths_at_vtimes
    depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m
    missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes
    x_velocities = ensure_numeric(x_velocities)
    x_velocities = x_velocities * -1.0 # Swap axis around
    y_velocities = ensure_numeric(y_velocities)
    y_velocities = y_velocities * -1.0 # Swap axis around
    
    fid = open(out_file,'w')

    assert len(depths_at_vtimes) == len(vtimes)
    start_time = None
    #start_time = 0.0
    #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes,
    #                                           velocities):
    for i in xrange(len(vtimes)):
        if not depths_at_vtimes[i] == missing:
            # Make the times start at zero.
            if start_time is None:
                start_time = vtimes[i]            
            fid.write(str(vtimes[i]-start_time) \
                      + ',' + str(depths_at_vtimes[i]) \
                      + ',' + str(x_velocities[i]) \
                      + ',' + str(y_velocities[i])+'\n')

    fid.close()
    
    # Since there is a new time reference save the depth info using this
    # new reference.
    fid = open(depth_file[:-4] + '_new_time'+depth_file[-4:],'w')
    sensors[:,0] -= start_time
    #print "depth_file", depth_file
    #print "start_time", start_time
    for i in xrange(len(dtimes)):
        fid.write(str(sensors[i,0]))
        for j in xrange(1,len(sensors[0])):
            fid.write(' ' + str(sensors[i,j]))
        fid.write('\n')
    fid.close()

    
    
#-------------------------------------------------------------------
if __name__ == "__main__":
    from os import getenv
    from os.path import join
    home = getenv('INUNDATIONHOME')
    Hinwood_dir = join(home,'data','flumes','Hinwood_2008')
    raw_data_dir = join(Hinwood_dir, 'raw_data')

    # Test 1 Run 5
    combine_velocity_depth(join(raw_data_dir,'T1R5velfilt.txt'),
                           join(raw_data_dir,'T1R5pressfilt.txt'),
                           join(Hinwood_dir, 'T1R5_boundary.csv'))    
    # Create the tsm file
    prepare_time_boundary(join(Hinwood_dir, 'T1R5_boundary.tsm'))
    
    # Test 2 Run 7
    combine_velocity_depth(join(raw_data_dir,'T2R7velfilt.txt'),
                           join(raw_data_dir,'T2R7pressfilt.txt'),
                           join(Hinwood_dir, 'T2R7_boundary.csv'))
    # Create the tsm file
    prepare_time_boundary(join(Hinwood_dir, 'T2R7_boundary.tsm'))