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

"""

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

from os.path import join

# 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')
    
def csv2tms(filename, offshore_bed_elevation):
    """Convert benchmark 2 time series to NetCDF tms file.
    the filename is the name of the output tms file, eg 'hi.tsm'.
    There must be an equivalent .csv file, eg 'hi.csv'.
   
    """



    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] = float(fields[1])
        depth = Q[i] - offshore_bed_elevation
        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()

def prepare_time_boundary(metadata_dic, raw_data_dir, output_dir):
    """
    """
    scenario_id = metadata_dic['scenario_id']
    velocity_file = join(raw_data_dir,scenario_id+'velfilt.txt')
    depth_file = join(raw_data_dir,scenario_id+'pressfilt.txt')
    out_file = join(output_dir, scenario_id+'_boundary.csv')
    
    combine_velocity_depth(velocity_file, depth_file, out_file)  
    tsm_file = out_file[:-4] + '.tsm'
    if metadata_dic['xleft'][1] >= 0.0:
        # This should be a -ve value, since the still water level is the
        # z origin.
        print "Warning: The z origin seems incorrect."
        
    csv2tms(tsm_file, metadata_dic['xleft'][1])
    
    
#-------------------------------------------------------------------
if __name__ == "__main__":
    pass