"""Create mesh and time boundary for the Okushiri island validation
"""


from Numeric import array, zeros, Float, allclose

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

import project_truescale

def prepare_bathymetry(filename):
    """Convert benchmark 2 bathymetry to NetCDF pts file.
    This is a 'throw-away' code taylor made for files like
    'okushiri_truescale_bathymetry.txt' from the LWRU2 benchmark
    """

    print 'Creating', filename
    
    # Read the ascii (.txt) version of this file,
    # make it comma separated and invert the bathymetry
    # (Below mean sea level should be negative)
    infile = open(filename[:-4] + '.txt')

    points = []
    attribute = []
    for line in infile.readlines()[1:]: #Skip first line (the header)
        fields = line.strip().split()

        x = float(fields[0])
        y = float(fields[1])
        z = -float(fields[2]) # Bathymetry is inverted in original file
        
        points.append([x,y])
        attribute.append(z)
    infile.close()

    # Convert to geospatial data and store as NetCDF
    G = Geospatial_data(data_points=points,
                        attributes=attribute)
    G.export_points_file(filename)
    


def prepare_timeboundary(filename):
    """Convert benchmark 2 time series to NetCDF tms file.
    This is a 'throw-away' code taylor made for files like
    'okushiri_truescale_input.txt' from the LWRU2 benchmark
    """

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


    print 'Creating', filename

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

    # Skip first line
    line = fid.readline()

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


    N = len(lines)
    T = zeros(N, Float)  #Time
    Q = zeros(N, Float)  #Values

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

        T[i] = float(fields[0])
        Q[i] = float(fields[1])


    # Create tms NetCDF file

    fid = NetCDFFile(filename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for truescale okushiri wavetank experiment'
    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'][:] = 0.0

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

    fid.close()


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


    # Prepare bathymetry
    prepare_bathymetry(project_truescale.bathymetry_filename)

    # Prepare time boundary
    prepare_timeboundary(project_truescale.boundary_filename)