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

"""
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 os.path import join
  
from anuga.utilities.numerical_tools import ensure_numeric

from interp import interp


import project

# 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, metadata_dic):
    """
    
    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, cm/sec],
    [y-velocity], [z-velocity]

    The format for the pressure file is
    [time, sec], [mm above SWL for sensor A], many other sensors...

    The format of the output file is;
    [time, sec], [m above SWL for sensor A],
    [x-velocity +ve is towards the shore, m/sec],
    [y-velocity +ve is towards the left of the tank, looking from the
    generator to the shore, m/sec]
    
    
    """
    missing = 1e+20
    
    # Read velocity file
    #print "*********************"
    #print "velocity_file", 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
    #print "depth_file", 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)
    #print "metadata_dic['scenario_id']", metadata_dic['scenario_id']
 #    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)
    # Swap axis around convert cm/sec to m/sec
    x_velocities = x_velocities * -0.01 
    y_velocities = ensure_numeric(y_velocities)
    # Swap axis around convert cm/sec to m/sec
    y_velocities = y_velocities * -0.01
    
    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]
            final_time = vtimes[i]-start_time
            fid.write(str(final_time) \
                      + ',' + str(depths_at_vtimes[i]) \
                      + ',' + str(x_velocities[i]) \
                      + ',' + str(y_velocities[i])+'\n')

    fid.close()
    print "The start time for is", start_time
    if not start_time == metadata_dic['ANUGA_start_time']:
        raise ValueError, "The calc'ed and recorded start times are different"
    # Modify the sensor array to reflect the new start time
    sensors[:,0] -= start_time
    
    # Since there is a new time reference save the depth info using this
    # new reference.
    fid = open(depth_file[:-4] + '_exp_depth.csv','w')
    #print "depth_file", depth_file
    #print "start_time", start_time

    
    bed_elevation_list = metadata_dic['gauge_bed_elevation']
    # +2 due to Time column and gauge A
    max_j = len(bed_elevation_list)+2
    # Write a header
    fid.write('Time')
    for gauge_x in metadata_dic['gauge_x']:
        fid.write(',' + str(gauge_x)) 
    fid.write('\n')
    for i in xrange(len(dtimes)):       
        fid.write(str(sensors[i,0])) # Time
        # Don't write sensor A.
        # It is the boundary condition.
        for j, bed_elevation in map(None,
                                    xrange(2,max_j),
                                    bed_elevation_list):
            # depth, m
            gauge = sensors[i,j]/1000 - bed_elevation # Convert from mm to m
            fid.write(',' + str(gauge)) 
        fid.write('\n')
    fid.close()

    
    # Since there is a new time reference save the stage info using this
    # new reference.
    fid = open(depth_file[:-4] + '_exp_stage.csv','w')
    
    # Write a header
    fid.write('Time')
    for gauge_x in metadata_dic['gauge_x']:
        fid.write(',' + str(gauge_x)) 
    fid.write('\n')
    for i in xrange(len(dtimes)):       
        fid.write(str(sensors[i,0])) # Time
        # Don't write sensor A.
        # It is the boundary condition.
        for j in xrange(2,max_j):
            # stage, m
            gauge = sensors[i,j]/1000 # Convert from mm to m
            fid.write(',' + str(gauge)) 
        fid.write('\n')
    fid.close()
    
    return final_time

def prepare_time_boundary(metadata_dic, raw_data_dir, output_dir):
    """
    Use this if a project instance has already been created.
    """
    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')
    #out_depth_file = join(output_dir, scenario_id+'_exp_depth.csv')
    
    final_time = combine_velocity_depth(velocity_file, depth_file, out_file,
                                        metadata_dic)
    #print "metadata_dic['xleft'][1]", metadata_dic['xleft'][1]
    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."     
    tsm_file = out_file[:-4] + '.tsm'  
    csv2tms(tsm_file, metadata_dic['xleft'][1])
    return final_time

# Don't do this, since run-dam changes the  metadata_dic['xleft'][1],
# which is used by this function
def prepare_time_boundary_for_scenarios_warning(scenarios,
                                        outputdir_tag):

    for run_data in scenarios:       
        id = run_data['scenario_id']
        outputdir_name = id + outputdir_tag
        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
                                       outputdir_name=outputdir_name)
        prepare_time_boundary(run_data,
                              pro_instance.raw_data_dir,
                              pro_instance.boundarydir)
        
    
#-------------------------------------------------------------------
if __name__ == "__main__":
    
    import sys; sys.exit()
    
    from scenarios import scenarios
    outputdir_tag = "XXXX"
    prepare_time_boundary_for_scenarios(scenarios,
                                        outputdir_tag)