source: trunk/anuga_work/development/Rudy_vandrie_2007/blackbutt/create_hydrograph.py @ 7924

"""Create mesh and time boundary for Hydrograph example
"""


from Numeric import array, zeros, Float, allclose

from anuga.pmesh.mesh import *
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.coordinate_transforms.geo_reference import Geo_reference
from anuga.geospatial_data import Geospatial_data

import project


def prepare_hydrograph(filename):
    """Convert .asc hydrograph file to NetCDF tms file.
    This is a 'throw-away' code taylor made for this type of file

    Two columns
    Time in hours    flow [m^3/s]
    """

    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import array

    assert filename[-4:] == '.tms'
   
    outfilename = filename
    infilename = filename[:-4] + '.asc'

    print 'Creating', outfilename

    # Read the ascii (.txt) version of this file
    fid = open(infilename)

    # Read all lines and search for selected hydrograph
    lines = fid.readlines()
    fid.close()

    time = []
    flow = []
    for line in lines:
        fields = line.split()
        time.append( float(fields[0])*3600   ) # Convert h to s
        flow.append( float(fields[1]) )        


    # Convert to NetCDF
    N = len(time)
    T = array(time, Float)  # Time (seconds)
    Q = array(flow, Float)  # Values (m^3/s)


    # Create tms NetCDF file
    fid = NetCDFFile(outfilename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Hydrograph example'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', Float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

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

    fid.close()


#-------------------------------------------------------------
if __name__ == "__main__":


    # Prepare hydrograph
    prepare_hydrograph(filename='Q/QPMF_Rot_Sub13.tms')