source: anuga_work/development/Rudy_vandrie_2007/create_hydro_example.py @ 6409

"""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, hydrograph, field):
    """Convert .out file to NetCDF tms file.
    This is a 'throw-away' code taylor made for this type of file

    hydrograph is an identifier e.g.
    HYDROGRAPHS_SUB17       

    field is one of the column headers.

    Time will always be read and converted into seconds.

    The file has contents like

    #####START_HYDROGRAPHS_SUB17       
     Time    Rain  Rainperv    Qtop      Qbot      Qper      Qimp    Qinto_OS   Qout_OS  Stage
      0.0    0.00    0.00     0.000     0.000     0.000     0.000     0.000     0.000     0.000
      5.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
     10.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
     15.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000

    
    """

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

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

    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()

    found = False
    search_string = 'START_'+hydrograph
    print 'looking for', search_string
    for i, line in enumerate(lines):
        if line.startswith('####'):
            if line.find(search_string) >= 0:
                print 'Found', line
                found = True
                break

    if found is False:
        msg = 'Did not find', hydrograph
        raise Exception, msg

    # Select data column
    headers = lines[i+1]
    fields = headers.split()
    index = fields.index(field)
    print 'Found header #%d: "%s"' %(index, fields[index])


    # Read data
    search_string = 'END_'+hydrograph
    time = []
    data = []    
    for line in lines[i+2:]:
        if line.startswith('####'):
            if line.find(search_string) >= 0:
                break        
        
        fields = line.split()
        time.append(float(fields[0]))
        data.append(float(fields[index]))


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

    # At this point, the hydrograph is stored in arrays T (time) and Q
    

    # 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=project.boundary_filename,
                       hydrograph='HYDROGRAPHS_SUB1',
                       field='Qinto_OS')