source: development/pyvolution-1d/data_manager.py @ 3322

"""datamanager.py - input output for AnuGA


This module takes care of reading and writing datafiles such as topograhies,
model output, etc


Formats used within AnuGA:

.sww: Netcdf format for storing model output f(t,x,y)
.tms: Netcdf format for storing time series f(t)

.xya: ASCII format for storing arbitrary points and associated attributes
.pts: NetCDF format for storing arbitrary points and associated attributes

.asc: ASCII format of regular DEMs as output from ArcView
.prj: Associated ArcView file giving more meta data for asc format
.ers: ERMapper header format of regular DEMs for ArcView

.dem: NetCDF representation of regular DEM data

.tsh: ASCII format for storing meshes and associated boundary and region info
.msh: NetCDF format for storing meshes and associated boundary and region info

.nc: Native ferret NetCDF format
.geo: Houdinis ascii geometry format (?)


A typical dataflow can be described as follows

Manually created files:
ASC, PRJ:     Digital elevation models (gridded)
TSH:          Triangular meshes (e.g. created from pmesh)
NC            Model outputs for use as boundary conditions (e.g from MOST)


AUTOMATICALLY CREATED FILES:

ASC, PRJ  ->  DEM  ->  PTS: Conversion of DEM's to native pts file

NC -> SWW: Conversion of MOST bundary files to boundary sww

PTS + TSH -> TSH with elevation: Least squares fit

TSH -> SWW:  Conversion of TSH to sww viewable using Swollen

TSH + Boundary SWW -> SWW: Simluation using pyvolution

"""
import os

from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
     searchsorted, zeros, allclose, around


from coordinate_transforms.geo_reference import Geo_reference

def make_filename(s):
    """Transform argument string into a suitable filename
    """

    s = s.strip()
    s = s.replace(' ', '_')
    s = s.replace('(', '')
    s = s.replace(')', '')
    s = s.replace('__', '_')

    return s


def check_dir(path, verbose=None):
    """Check that specified path exists.
    If path does not exist it will be created if possible

    USAGE:
       checkdir(path, verbose):

    ARGUMENTS:
        path -- Directory
        verbose -- Flag verbose output (default: None)

    RETURN VALUE:
        Verified path including trailing separator

    """

    import os, sys
    import os.path

    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
        unix = 0
    else:
        unix = 1


    if path[-1] != os.sep:
        path = path + os.sep  # Add separator for directories

    path = os.path.expanduser(path) # Expand ~ or ~user in pathname
    if not (os.access(path,os.R_OK and os.W_OK) or path == ''):
        try:
            exitcode=os.mkdir(path)

            # Change access rights if possible
            #
            if unix:
                exitcode=os.system('chmod 775 '+path)
            else:
                pass  # FIXME: What about acces rights under Windows?

            if verbose: print 'MESSAGE: Directory', path, 'created.'

        except:
            print 'WARNING: Directory', path, 'could not be created.'
            if unix:
                path = '/tmp/'
            else:
                path = 'C:'

            print 'Using directory %s instead' %path

    return(path)



def del_dir(path):
    """Recursively delete directory path and all its contents
    """

    import os

    if os.path.isdir(path):
        for file in os.listdir(path):
            X = os.path.join(path, file)


            if os.path.isdir(X) and not os.path.islink(X):
                del_dir(X)
            else:
                try:
                    os.remove(X)
                except:
                    print "Could not remove file %s" %X

        os.rmdir(path)



def create_filename(datadir, filename, format, size=None, time=None):

    import os
    #from config import data_dir

    FN = check_dir(datadir) + filename

    if size is not None:
        FN += '_size%d' %size

    if time is not None:
        FN += '_time%.2f' %time

    FN += '.' + format
    return FN


def get_files(datadir, filename, format, size):
    """Get all file (names) with gven name, size and format
    """

    import glob

    import os
    #from config import data_dir

    dir = check_dir(datadir)

    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
    return glob.glob(pattern)



#Generic class for storing output to e.g. visualisation or checkpointing
class Data_format:
    """Generic interface to data formats
    """


    def __init__(self, domain, extension, mode = 'w'):
        assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\
                                        '''   'w' (write)'''+\
                                        '''   'r' (read)''' +\
                                        '''   'a' (append)'''

        #Create filename
        #self.filename = create_filename(domain.get_datadir(),
        #domain.get_name(), extension, len(domain))


        self.filename = create_filename(domain.get_datadir(),
                        domain.get_name(), extension)

        #print 'F', self.filename
        self.timestep = 0
        #I CHANGED THIS
        #self.number_of_volumes = len(domain)
        self.number_of_volumes = domain.number_of_elements
        self.domain = domain


        #FIXME: Should we have a general set_precision function?



#Class for storing output to e.g. visualisation
class Data_format_sww(Data_format):
    """Interface to native NetCDF format (.sww) for storing model output

    There are two kinds of data

    1: Constant data: Vertex coordinates and field values. Stored once
    2: Variable data: Conserved quantities. Stored once per timestep.

    All data is assumed to reside at vertex locations.
    """


    def __init__(self, domain, mode = 'w',\
                 max_size = 2000000000,
                 recursion = False):
        from Scientific.IO.NetCDF import NetCDFFile
        from Numeric import Int, Float, Float32

        self.precision = Float32 #Use single precision
        if hasattr(domain, 'max_size'):
            self.max_size = domain.max_size #file size max is 2Gig
        else:
            self.max_size = max_size
        self.recursion = recursion
        self.mode = mode

        Data_format.__init__(self, domain, 'sww', mode)


        # NetCDF file definition
        fid = NetCDFFile(self.filename, mode)

        if mode == 'w':

            #Create new file
            fid.institution = 'Geoscience Australia'
            fid.description = 'Output from pyvolution suitable for plotting'

            if domain.smooth:
                fid.smoothing = 'Yes'
            else:
                fid.smoothing = 'No'

            fid.order = domain.default_order

            if hasattr(domain, 'texture'):
                fid.texture = domain.texture
            #else:
            #    fid.texture = 'None'

            #Reference point
            #Start time in seconds since the epoch (midnight 1/1/1970)
            #FIXME: Use Georef
            fid.starttime = domain.starttime

            # dimension definitions
            fid.createDimension('number_of_volumes', self.number_of_volumes)
            #fid.createDimension('number_of_vertices', 3)
            fid.createDimension('number_of_vertices', 2)

            if domain.smooth is True:
                #i changed this
                #fid.createDimension('number_of_points', len(domain.vertexlist))
                fid.createDimension('number_of_points', domain.number_of_elements+1)
            else:
                #fid.createDimension('number_of_points', 3*self.number_of_volumes)
                fid.createDimension('number_of_points', 2*self.number_of_volumes)

            fid.createDimension('number_of_timesteps', None) #extensible

            # variable definitions
            fid.createVariable('x', self.precision, ('number_of_points',))
            #fid.createVariable('y', self.precision, ('number_of_points',))
            fid.createVariable('elevation', self.precision, ('number_of_points',))
            if domain.geo_reference is not None:
                domain.geo_reference.write_NetCDF(fid)

            #FIXME: Backwards compatibility
            fid.createVariable('z', self.precision, ('number_of_points',))
            #################################

            fid.createVariable('volumes', Int, ('number_of_volumes',
                                                'number_of_vertices'))

            fid.createVariable('time', self.precision,
                               ('number_of_timesteps',))

            fid.createVariable('stage', self.precision,
                               ('number_of_timesteps',
                                'number_of_points'))

            fid.createVariable('xmomentum', self.precision,
                               ('number_of_timesteps',
                                'number_of_points'))

            #fid.createVariable('ymomentum', self.precision,
            #                   ('number_of_timesteps',
            #                    'number_of_points'))

        #Close
        fid.close()


    def store_connectivity(self):
        """Specialisation of store_connectivity for net CDF format

        Writes x,y,z coordinates of triangles constituting
        the bed elevation.
        """

        from Scientific.IO.NetCDF import NetCDFFile

        from Numeric import concatenate, Int

        domain = self.domain

        #Get NetCDF
        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append

        # Get the variables
        x = fid.variables['x']
        #y = fid.variables['y']
        z = fid.variables['elevation']

        volumes = fid.variables['volumes']

        # Get X, Y and bed elevation Z
        Q = domain.quantities['elevation']
        #X,Y,Z,V = Q.get_vertex_values(xy=True,
        #                             precision = self.precision)
        X,Z,V = Q.get_vertex_values(x=True,
                                      precision = self.precision)

        x[:] = X.astype(self.precision)
        #y[:] = Y.astype(self.precision)
        z[:] = Z.astype(self.precision)

        #FIXME: Backwards compatibility
        z = fid.variables['z']
        z[:] = Z.astype(self.precision)
        ################################

        volumes[:] = V.astype(volumes.typecode())

        #Close
        fid.close()

    def store_timestep(self, names):
        """Store time and named quantities to file
        """
        from Scientific.IO.NetCDF import NetCDFFile
        import types
        from time import sleep
        from os import stat

        #minimum_allowed_depth = 0.001
        minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
        from Numeric import choose

        #Get NetCDF
        retries = 0
        file_open = False
        while not file_open and retries < 10:
            try:
                fid = NetCDFFile(self.filename, 'a') #Open existing file
            except IOError:
                #This could happen if someone was reading the file.
                #In that case, wait a while and try again
                msg = 'Warning (store_timestep): File %s could not be opened'\
                      %self.filename
                msg += ' - trying step %s again' %self.domain.time
                print msg
                retries += 1
                sleep(1)
            else:
                file_open = True

        if not file_open:
            msg = 'File %s could not be opened for append' %self.filename
            raise msg



        #Check to see if the file is already too big:
        time = fid.variables['time']
        i = len(time)+1
        file_size = stat(self.filename)[6]
        file_size_increase =  file_size/i
        if file_size + file_size_increase > self.max_size*(2**self.recursion):
            #in order to get the file name and start time correct,
            #I change the domian.filename and domain.starttime.
            #This is the only way to do this without changing
            #other modules (I think).

            #write a filename addon that won't break swollens reader
            #(10.sww is bad)
            filename_ext = '_time_%s'%self.domain.time
            filename_ext = filename_ext.replace('.', '_')
            #remember the old filename, then give domain a
            #name with the extension
            old_domain_filename = self.domain.filename
            if not self.recursion:
                self.domain.filename = self.domain.filename+filename_ext

            #change the domain starttime to the current time
            old_domain_starttime = self.domain.starttime
            self.domain.starttime = self.domain.time

            #build a new data_structure.
            next_data_structure=\
                Data_format_sww(self.domain, mode=self.mode,\
                                max_size = self.max_size,\
                                recursion = self.recursion+1)
            if not self.recursion:
                print '    file_size = %s'%file_size
                print '    saving file to %s'%next_data_structure.filename
            #set up the new data_structure
            self.domain.writer = next_data_structure

            #FIXME - could be cleaner to use domain.store_timestep etc.
            next_data_structure.store_connectivity()
            next_data_structure.store_timestep(names)
            fid.sync()
            fid.close()

            #restore the old starttime and filename
            self.domain.starttime = old_domain_starttime
            self.domain.filename = old_domain_filename
        else:
            self.recursion = False
            domain = self.domain

            # Get the variables
            time = fid.variables['time']
            stage = fid.variables['stage']
            xmomentum = fid.variables['xmomentum']
            #ymomentum = fid.variables['ymomentum']
            i = len(time)

            #Store time
            time[i] = self.domain.time


            if type(names) not in [types.ListType, types.TupleType]:
                names = [names]

            for name in names:
                # Get quantity
                Q = domain.quantities[name]
                #A,V = Q.get_vertex_values(xy = False,
                #                          precision = self.precision)
                A,V = Q.get_vertex_values(x = False,
                                          precision = self.precision)

                #FIXME: Make this general (see below)
                if name == 'stage':
                    z = fid.variables['elevation']
                    #print z[:]
                    #print A-z[:]
                    A = choose( A-z[:] >= minimum_allowed_depth, (z[:], A))
                    stage[i,:] = A.astype(self.precision)
                elif name == 'xmomentum':
                    xmomentum[i,:] = A.astype(self.precision)
                elif name == 'ymomentum':
                    ymomentum[i,:] = A.astype(self.precision)

                #As in....
                #eval( name + '[i,:] = A.astype(self.precision)' )
                #FIXME: But we need a UNIT test for that before refactoring



            #Flush and close
            fid.sync()
            fid.close()



#Class for handling checkpoints data
class Data_format_cpt(Data_format):
    """Interface to native NetCDF format (.cpt)
    """


    def __init__(self, domain, mode = 'w'):
        from Scientific.IO.NetCDF import NetCDFFile
        from Numeric import Int, Float, Float

        self.precision = Float #Use full precision

        Data_format.__init__(self, domain, 'sww', mode)


        # NetCDF file definition
        fid = NetCDFFile(self.filename, mode)

        if mode == 'w':
            #Create new file
            fid.institution = 'Geoscience Australia'
            fid.description = 'Checkpoint data'
            #fid.smooth = domain.smooth
            fid.order = domain.default_order

            # dimension definitions
            fid.createDimension('number_of_volumes', self.number_of_volumes)
            fid.createDimension('number_of_vertices', 3)

            #Store info at all vertices (no smoothing)
            fid.createDimension('number_of_points', 3*self.number_of_volumes)
            fid.createDimension('number_of_timesteps', None) #extensible

            # variable definitions

            #Mesh
            fid.createVariable('x', self.precision, ('number_of_points',))
            fid.createVariable('y', self.precision, ('number_of_points',))


            fid.createVariable('volumes', Int, ('number_of_volumes',
                                                'number_of_vertices'))

            fid.createVariable('time', self.precision,
                               ('number_of_timesteps',))

            #Allocate space for all quantities
            for name in domain.quantities.keys():
                fid.createVariable(name, self.precision,
                                   ('number_of_timesteps',
                                    'number_of_points'))

        #Close
        fid.close()


    def store_checkpoint(self):
        """
        Write x,y coordinates of triangles.
        Write connectivity (
        constituting
        the bed elevation.
        """

        from Scientific.IO.NetCDF import NetCDFFile

        from Numeric import concatenate

        domain = self.domain

        #Get NetCDF
        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append

        # Get the variables
        x = fid.variables['x']
        y = fid.variables['y']

        volumes = fid.variables['volumes']

        # Get X, Y and bed elevation Z
        Q = domain.quantities['elevation']
        X,Y,Z,V = Q.get_vertex_values(xy=True,
                                      precision = self.precision)



        x[:] = X.astype(self.precision)
        y[:] = Y.astype(self.precision)
        z[:] = Z.astype(self.precision)

        volumes[:] = V

        #Close
        fid.close()


    def store_timestep(self, name):
        """Store time and named quantity to file
        """
        from Scientific.IO.NetCDF import NetCDFFile
        from time import sleep

        #Get NetCDF
        retries = 0
        file_open = False
        while not file_open and retries < 10:
            try:
                fid = NetCDFFile(self.filename, 'a') #Open existing file
            except IOError:
                #This could happen if someone was reading the file.
                #In that case, wait a while and try again
                msg = 'Warning (store_timestep): File %s could not be opened'\
                          %self.filename
                msg += ' - trying again'
                print msg
                retries += 1
                sleep(1)
            else:
                file_open = True

            if not file_open:
                msg = 'File %s could not be opened for append' %self.filename
            raise msg


        domain = self.domain

        # Get the variables
        time = fid.variables['time']
        stage = fid.variables['stage']
        i = len(time)

        #Store stage
        time[i] = self.domain.time

        # Get quantity
        Q = domain.quantities[name]
        A,V = Q.get_vertex_values(xy=False,
                                  precision = self.precision)

        stage[i,:] = A.astype(self.precision)

        #Flush and close
        fid.sync()
        fid.close()





#Function for storing xya output
#FIXME Not done yet for this version
class Data_format_xya(Data_format):
    """Generic interface to data formats
    """


    def __init__(self, domain, mode = 'w'):
        from Scientific.IO.NetCDF import NetCDFFile
        from Numeric import Int, Float, Float32

        self.precision = Float32 #Use single precision

        Data_format.__init__(self, domain, 'xya', mode)



    #FIXME -This is the old xya format
    def store_all(self):
        """Specialisation of store all for xya format

        Writes x,y,z coordinates of triangles constituting
        the bed elevation.
        """

        from Numeric import concatenate

        domain = self.domain

        fd = open(self.filename, 'w')


        if domain.smooth is True:
            number_of_points =  len(domain.vertexlist)
        else:
            number_of_points = 3*self.number_of_volumes

        numVertAttrib = 3 #Three attributes is what is assumed by the xya format

        fd.write(str(number_of_points) + " " + str(numVertAttrib) +\
                 " # <vertex #> <x> <y> [attributes]" + "\n")


        # Get X, Y, bed elevation and friction (index=0,1)
        X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values',
                                           indices = (0,1), precision = self.precision)

        bed_eles = A[:,0]
        fricts = A[:,1]

        # Get stage (index=0)
        B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities',
                                       indices = (0,), precision = self.precision)

        stages = B[:,0]

        #<vertex #> <x> <y> [attributes]
        for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles,
                                               stages, fricts):

            s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict)
            fd.write(s)

        #close
        fd.close()


    def store_timestep(self, t, V0, V1, V2):
        """Store time, water heights (and momentums) to file
        """
        pass


#Auxiliary
def write_obj(filename,x,y,z):
    """Store x,y,z vectors into filename (obj format)
       Vectors are assumed to have dimension (M,3) where
       M corresponds to the number elements.
       triangles are assumed to be disconnected

       The three numbers in each vector correspond to three vertices,

       e.g. the x coordinate of vertex 1 of element i is in x[i,1]

    """
    #print 'Writing obj to %s' % filename

    import os.path

    root, ext = os.path.splitext(filename)
    if ext == '.obj':
        FN = filename
    else:
        FN = filename + '.obj'


    outfile = open(FN, 'wb')
    outfile.write("# Triangulation as an obj file\n")

    M, N = x.shape
    assert N==3  #Assuming three vertices per element

    for i in range(M):
        for j in range(N):
            outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j]))

    for i in range(M):
        base = i*N
        outfile.write("f %d %d %d\n" % (base+1,base+2,base+3))

    outfile.close()


#########################################################
#Conversion routines
########################################################

def sww2obj(basefilename, size):
    """Convert netcdf based data output to obj
    """
    from Scientific.IO.NetCDF import NetCDFFile

    from Numeric import Float, zeros

    #Get NetCDF
    FN = create_filename('.', basefilename, 'sww', size)
    print 'Reading from ', FN
    fid = NetCDFFile(FN, 'r')  #Open existing file for read


    # Get the variables
    x = fid.variables['x']
    y = fid.variables['y']
    z = fid.variables['elevation']
    time = fid.variables['time']
    stage = fid.variables['stage']

    M = size  #Number of lines
    xx = zeros((M,3), Float)
    yy = zeros((M,3), Float)
    zz = zeros((M,3), Float)

    for i in range(M):
        for j in range(3):
            xx[i,j] = x[i+j*M]
            yy[i,j] = y[i+j*M]
            zz[i,j] = z[i+j*M]

    #Write obj for bathymetry
    FN = create_filename('.', basefilename, 'obj', size)
    write_obj(FN,xx,yy,zz)


    #Now read all the data with variable information, combine with
    #x,y info and store as obj

    for k in range(len(time)):
        t = time[k]
        print 'Processing timestep %f' %t

        for i in range(M):
            for j in range(3):
                zz[i,j] = stage[k,i+j*M]


        #Write obj for variable data
        #FN = create_filename(basefilename, 'obj', size, time=t)
        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
        write_obj(FN,xx,yy,zz)


def dat2obj(basefilename):
    """Convert line based data output to obj
    FIXME: Obsolete?
    """

    import glob, os
    from config import data_dir


    #Get bathymetry and x,y's
    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()

    from Numeric import zeros, Float

    M = len(lines)  #Number of lines
    x = zeros((M,3), Float)
    y = zeros((M,3), Float)
    z = zeros((M,3), Float)

    ##i = 0
    for i, line in enumerate(lines):
        tokens = line.split()
        values = map(float,tokens)

        for j in range(3):
            x[i,j] = values[j*3]
            y[i,j] = values[j*3+1]
            z[i,j] = values[j*3+2]

        ##i += 1


    #Write obj for bathymetry
    write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z)


    #Now read all the data files with variable information, combine with
    #x,y info
    #and store as obj

    files = glob.glob(data_dir+os.sep+basefilename+'*.dat')

    for filename in files:
        print 'Processing %s' % filename

        lines = open(data_dir+os.sep+filename,'r').readlines()
        assert len(lines) == M
        root, ext = os.path.splitext(filename)

        #Get time from filename
        i0 = filename.find('_time=')
        if i0 == -1:
            #Skip bathymetry file
            continue

        i0 += 6  #Position where time starts
        i1 = filename.find('.dat')

        if i1 > i0:
            t = float(filename[i0:i1])
        else:
            raise 'Hmmmm'



        ##i = 0
        for i, line in enumerate(lines):
            tokens = line.split()
            values = map(float,tokens)

            for j in range(3):
                z[i,j] = values[j]

            ##i += 1

        #Write obj for variable data
        write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z)


def filter_netcdf(filename1, filename2, first=0, last=None, step = 1):
    """Read netcdf filename1, pick timesteps first:step:last and save to
    nettcdf file filename2
    """
    from Scientific.IO.NetCDF import NetCDFFile

    #Get NetCDF
    infile = NetCDFFile(filename1, 'r')  #Open existing file for read
    outfile = NetCDFFile(filename2, 'w')  #Open new file


    #Copy dimensions
    for d in infile.dimensions:
        outfile.createDimension(d, infile.dimensions[d])

    for name in infile.variables:
        var = infile.variables[name]
        outfile.createVariable(name, var.typecode(), var.dimensions)


    #Copy the static variables
    for name in infile.variables:
        if name == 'time' or name == 'stage':
            pass
        else:
            #Copy
            outfile.variables[name][:] = infile.variables[name][:]

    #Copy selected timesteps
    time = infile.variables['time']
    stage = infile.variables['stage']

    newtime = outfile.variables['time']
    newstage = outfile.variables['stage']

    if last is None:
        last = len(time)

    selection = range(first, last, step)
    for i, j in enumerate(selection):
        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
        newtime[i] = time[j]
        newstage[i,:] = stage[j,:]

    #Close
    infile.close()
    outfile.close()


#Get data objects
def get_dataobject(domain, mode='w'):
    """Return instance of class of given format using filename
    """

    cls = eval('Data_format_%s' %domain.format)
    return cls(domain, mode)

def xya2pts(basename_in, basename_out=None, verbose=False,
            #easting_min=None, easting_max=None,
            #northing_min=None, northing_max=None,
            stride = 1,
            attribute_name = 'elevation',
            z_func = None):
    """Read points data from ascii (.xya)

    Example:

              x(m)        y(m)        z(m)
         0.00000e+00  0.00000e+00  1.3535000e-01
         0.00000e+00  1.40000e-02  1.3535000e-01



    Convert to NetCDF pts format which is

    points:  (Nx2) Float array
    elevation: N Float array

    Only lines that contain three numeric values are processed

    If z_func is specified, it will be applied to the third field
    """

    import os
    #from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, arrayrange, concatenate

    root, ext = os.path.splitext(basename_in)

    if ext == '': ext = '.xya'

    #Get NetCDF
    infile = open(root + ext, 'r')  #Open existing xya file for read

    if verbose: print 'Reading xya points from %s' %(root + ext)

    points = []
    attribute = []
    for i, line in enumerate(infile.readlines()):

        if i % stride != 0: continue

        fields = line.split()

        try:
            assert len(fields) == 3
        except:
            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line)

        try:
            x = float( fields[0] )
            y = float( fields[1] )
            z = float( fields[2] )
        except:
            continue

        points.append( [x, y] )

        if callable(z_func):
            attribute.append(z_func(z))
        else:
            attribute.append(z)


    #Get output file
    if basename_out == None:
        ptsname = root + '.pts'
    else:
        ptsname = basename_out + '.pts'

    if verbose: print 'Store to NetCDF file %s' %ptsname
    write_ptsfile(ptsname, points, attribute, attribute_name)



######Obsoleted by export_points in load_mesh
def write_ptsfile(ptsname, points, attribute, attribute_name = None,
                  zone=None, xllcorner=None, yllcorner=None):
    """Write points and associated attribute to pts (NetCDF) format
    """

    print 'WARNING: write_ptsfile is obsolete. Use export_points from load_mesh.loadASCII instead.'

    from Numeric import Float

    if attribute_name is None:
        attribute_name = 'attribute'


    from Scientific.IO.NetCDF import NetCDFFile

    # NetCDF file definition
    outfile = NetCDFFile(ptsname, 'w')


    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF pts format for compact and '\
                          'portable storage of spatial point data'


    #Georeferencing
    from coordinate_transforms.geo_reference import Geo_reference
    if zone is None:
        assert xllcorner is None, 'xllcorner must be None'
        assert yllcorner is None, 'yllcorner must be None'        
        Geo_reference().write_NetCDF(outfile)
    else:
        Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile)        
        


    outfile.createDimension('number_of_points', len(points))
    outfile.createDimension('number_of_dimensions', 2) #This is 2d data

    # variable definitions
    outfile.createVariable('points', Float, ('number_of_points',
                                             'number_of_dimensions'))
    outfile.createVariable(attribute_name, Float, ('number_of_points',))

    # Get handles to the variables
    nc_points = outfile.variables['points']
    nc_attribute = outfile.variables[attribute_name]

    #Store data
    nc_points[:, :] = points
    nc_attribute[:] = attribute

    outfile.close()


def dem2pts(basename_in, basename_out=None, verbose=False,
            easting_min=None, easting_max=None,
            northing_min=None, northing_max=None):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Convert to NetCDF pts format which is

    points:  (Nx2) Float array
    elevation: N Float array
    """

    #FIXME: Can this be written feasibly using write_pts?

    import os
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, zeros, reshape

    root = basename_in

    #Get NetCDF
    infile = NetCDFFile(root + '.dem', 'r')  #Open existing netcdf file for read

    if verbose: print 'Reading DEM from %s' %(root + '.dem')

    ncols = infile.ncols[0]
    nrows = infile.nrows[0]
    xllcorner = infile.xllcorner[0]  #Easting of lower left corner
    yllcorner = infile.yllcorner[0]  #Northing of lower left corner
    cellsize = infile.cellsize[0]
    NODATA_value = infile.NODATA_value[0]
    dem_elevation = infile.variables['elevation']

    zone = infile.zone[0]
    false_easting = infile.false_easting[0]
    false_northing = infile.false_northing[0]

    #Text strings
    projection = infile.projection
    datum = infile.datum
    units = infile.units


    #Get output file
    if basename_out == None:
        ptsname = root + '.pts'
    else:
        ptsname = basename_out + '.pts'

    if verbose: print 'Store to NetCDF file %s' %ptsname
    # NetCDF file definition
    outfile = NetCDFFile(ptsname, 'w')

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
                      'of spatial point data'
    #assign default values
    if easting_min is None: easting_min = xllcorner
    if easting_max is None: easting_max = xllcorner + ncols*cellsize
    if northing_min is None: northing_min = yllcorner
    if northing_max is None: northing_max = yllcorner + nrows*cellsize

    #compute offsets to update georeferencing
    easting_offset = xllcorner - easting_min
    northing_offset = yllcorner - northing_min

    #Georeferencing
    outfile.zone = zone
    outfile.xllcorner = easting_min #Easting of lower left corner
    outfile.yllcorner = northing_min #Northing of lower left corner
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units


    #Grid info (FIXME: probably not going to be used, but heck)
    outfile.ncols = ncols
    outfile.nrows = nrows


    # dimension definitions
    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
    outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box)
    outfile.createDimension('number_of_dimensions', 2) #This is 2d data

    # variable definitions
    outfile.createVariable('points', Float, ('number_of_points',
                                             'number_of_dimensions'))
    outfile.createVariable('elevation', Float, ('number_of_points',))

    # Get handles to the variables
    points = outfile.variables['points']
    elevation = outfile.variables['elevation']

    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))

    #Store data
    #FIXME: Could perhaps be faster using array operations (Fixed 27/7/05)
    global_index = 0
    for i in range(nrows):
        if verbose and i%((nrows+10)/10)==0:
            print 'Processing row %d of %d' %(i, nrows)

        lower_index = global_index
        tpoints = zeros((ncols_in_bounding_box, 2), Float)
        telev =  zeros(ncols_in_bounding_box, Float)
        local_index = 0

        y = (nrows-i)*cellsize + yllcorner
        for j in range(ncols):

            x = j*cellsize + xllcorner
            if easting_min <= x <= easting_max and \
               northing_min <= y <= northing_max:
                tpoints[local_index, :] = [x-easting_min,y-northing_min]
                telev[local_index] = dem_elevation_r[i, j]
                global_index += 1
                local_index += 1

        upper_index = global_index

        if upper_index == lower_index + ncols_in_bounding_box:
            points[lower_index:upper_index, :] = tpoints
            elevation[lower_index:upper_index] = telev

    assert global_index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points'

    infile.close()
    outfile.close()



def _read_hecras_cross_sections(lines):
    """Return block of surface lines for each cross section
    Starts with SURFACE LINE,
    Ends with END CROSS-SECTION
    """

    points = []

    reading_surface = False
    for i, line in enumerate(lines):

        if len(line.strip()) == 0:    #Ignore blanks
            continue

        if lines[i].strip().startswith('SURFACE LINE'):
            reading_surface = True
            continue

        if lines[i].strip().startswith('END') and reading_surface:
            yield points
            reading_surface = False
            points = []

        if reading_surface:
            fields = line.strip().split(',')
            easting = float(fields[0])
            northing = float(fields[1])
            elevation = float(fields[2])
            points.append([easting, northing, elevation])




def hecras_cross_sections2pts(basename_in,
                              basename_out=None,
                              verbose=False):
    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)

    Example:


# RAS export file created on Mon 15Aug2005 11:42
# by HEC-RAS Version 3.1.1

BEGIN HEADER:
  UNITS: METRIC
  DTM TYPE: TIN
  DTM: v:\1\cit\perth_topo\river_tin
  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
  MAP PROJECTION: UTM
  PROJECTION ZONE: 50
  DATUM: AGD66
  VERTICAL DATUM:
  NUMBER OF REACHES:  19
  NUMBER OF CROSS-SECTIONS:  14206
END HEADER:


Only the SURFACE LINE data of the following form will be utilised

  CROSS-SECTION:
    STREAM ID:Southern-Wungong
    REACH ID:Southern-Wungong
    STATION:19040.*
    CUT LINE:
      405548.671603161 , 6438142.7594925
      405734.536092045 , 6438326.10404912
      405745.130459356 , 6438331.48627354
      405813.89633823 , 6438368.6272789
    SURFACE LINE:
     405548.67,   6438142.76,   35.37
     405552.24,   6438146.28,   35.41
     405554.78,   6438148.78,   35.44
     405555.80,   6438149.79,   35.44
     405559.37,   6438153.31,   35.45
     405560.88,   6438154.81,   35.44
     405562.93,   6438156.83,   35.42
     405566.50,   6438160.35,   35.38
     405566.99,   6438160.83,   35.37
     ...
   END CROSS-SECTION

    Convert to NetCDF pts format which is

    points:  (Nx2) Float array
    elevation: N Float array
    """

    #FIXME: Can this be written feasibly using write_pts?

    import os
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, zeros, reshape

    root = basename_in

    #Get ASCII file
    infile = open(root + '.sdf', 'r')  #Open SDF file for read

    if verbose: print 'Reading DEM from %s' %(root + '.sdf')

    lines = infile.readlines()
    infile.close()

    if verbose: print 'Converting to pts format'

    i = 0
    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
        i += 1

    assert lines[i].strip().upper() == 'BEGIN HEADER:'
    i += 1

    assert lines[i].strip().upper().startswith('UNITS:')
    units = lines[i].strip().split()[1]
    i += 1

    assert lines[i].strip().upper().startswith('DTM TYPE:')
    i += 1

    assert lines[i].strip().upper().startswith('DTM:')
    i += 1

    assert lines[i].strip().upper().startswith('STREAM')
    i += 1

    assert lines[i].strip().upper().startswith('CROSS')
    i += 1

    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
    projection = lines[i].strip().split(':')[1]
    i += 1

    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
    zone = int(lines[i].strip().split(':')[1])
    i += 1

    assert lines[i].strip().upper().startswith('DATUM:')
    datum = lines[i].strip().split(':')[1]
    i += 1

    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
    i += 1

    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
    i += 1

    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
    number_of_cross_sections = int(lines[i].strip().split(':')[1])
    i += 1


    #Now read all points
    points = []
    elevation = []
    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
        for k, entry in enumerate(entries):
            points.append(entry[:2])
            elevation.append(entry[2])


    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
          %(j+1, number_of_cross_sections)
    assert j+1 == number_of_cross_sections, msg

    #Get output file
    if basename_out == None:
        ptsname = root + '.pts'
    else:
        ptsname = basename_out + '.pts'

    #FIXME (DSG-ON): use loadASCII export_points_file
    if verbose: print 'Store to NetCDF file %s' %ptsname
    # NetCDF file definition
    outfile = NetCDFFile(ptsname, 'w')

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF pts format for compact and portable ' +\
                          'storage of spatial point data derived from HEC-RAS'

    #Georeferencing
    outfile.zone = zone
    outfile.xllcorner = 0.0
    outfile.yllcorner = 0.0
    outfile.false_easting = 500000    #FIXME: Use defaults from e.g. config.py
    outfile.false_northing = 1000000

    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units


    outfile.createDimension('number_of_points', len(points))
    outfile.createDimension('number_of_dimensions', 2) #This is 2d data

    # variable definitions
    outfile.createVariable('points', Float, ('number_of_points',
                                             'number_of_dimensions'))
    outfile.createVariable('elevation', Float, ('number_of_points',))

    # Get handles to the variables
    outfile.variables['points'][:] = points
    outfile.variables['elevation'][:] = elevation

    outfile.close()



def sww2dem(basename_in, basename_out = None,
            quantity = None,
            timestep = None,
            reduction = None,
            cellsize = 10,
            NODATA_value = -9999,
            easting_min = None,
            easting_max = None,
            northing_min = None,
            northing_max = None,
            expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned)
            verbose = False,
            origin = None,
            datum = 'WGS84',
            format = 'ers'):

    """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)

    Example (ASC):

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Also write accompanying file with same basename_in but extension .prj
    used to fix the UTM zone, datum, false northings and eastings.

    The prj format is assumed to be as

    Projection    UTM
    Zone          56
    Datum         WGS84
    Zunits        NO
    Units         METERS
    Spheroid      WGS84
    Xshift        0.0000000000
    Yshift        10000000.0000000000
    Parameters


    The parameter quantity must be the name of an existing quantity or
    an expression involving existing quantities. The default is
    'elevation'.

    if timestep (an index) is given, output quantity at that timestep

    if reduction is given use that to reduce quantity over all timesteps.

    datum

    format can be either 'asc' or 'ers'
    """

    import sys
    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
    from Numeric import array2string

    from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
    from util import apply_expression_to_dictionary

    msg = 'Format must be either asc or ers'
    assert format.lower() in ['asc', 'ers'], msg


    false_easting = 500000
    false_northing = 10000000

    if quantity is None:
        quantity = 'elevation'

    if reduction is None:
        reduction = max

    if basename_out is None:
        basename_out = basename_in + '_%s' %quantity

    swwfile = basename_in + '.sww'
    demfile = basename_out + '.' + format
    # Note the use of a .ers extension is optional (write_ermapper_grid will
    # deal with either option

    #Read sww file
    if verbose: print 'Reading from %s' %swwfile
    from Scientific.IO.NetCDF import NetCDFFile
    fid = NetCDFFile(swwfile)

    #Get extent and reference
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    volumes = fid.variables['volumes'][:]

    number_of_timesteps = fid.dimensions['number_of_timesteps']
    number_of_points = fid.dimensions['number_of_points']
    if origin is None:

        #Get geo_reference
        #sww files don't have to have a geo_ref
        try:
            geo_reference = Geo_reference(NetCDFObject=fid)
        except AttributeError, e:
            geo_reference = Geo_reference() #Default georef object

        xllcorner = geo_reference.get_xllcorner()
        yllcorner = geo_reference.get_yllcorner()
        zone = geo_reference.get_zone()
    else:
        zone = origin[0]
        xllcorner = origin[1]
        yllcorner = origin[2]



    #FIXME: Refactor using code from file_function.statistics
    #Something like print swwstats(swwname)
    if verbose:
        x = fid.variables['x'][:]
        y = fid.variables['y'][:]
        times = fid.variables['time'][:]
        print '------------------------------------------------'
        print 'Statistics of SWW file:'
        print '  Name: %s' %swwfile
        print '  Reference:'
        print '    Lower left corner: [%f, %f]'\
              %(xllcorner, yllcorner)
        print '    Start time: %f' %fid.starttime[0]
        print '  Extent:'
        print '    x [m] in [%f, %f], len(x) == %d'\
              %(min(x.flat), max(x.flat), len(x.flat))
        print '    y [m] in [%f, %f], len(y) == %d'\
              %(min(y.flat), max(y.flat), len(y.flat))
        print '    t [s] in [%f, %f], len(t) == %d'\
              %(min(times), max(times), len(times))
        print '  Quantities [SI units]:'
        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
            q = fid.variables[name][:].flat
            print '    %s in [%f, %f]' %(name, min(q), max(q))





    #Get quantity and reduce if applicable
    if verbose: print 'Processing quantity %s' %quantity

    #Turn NetCDF objects into Numeric arrays
    quantity_dict = {}
    for name in fid.variables.keys():
        quantity_dict[name] = fid.variables[name][:]


    q = apply_expression_to_dictionary(quantity, quantity_dict)



    if len(q.shape) == 2:
        #q has a time component and needs to be reduced along
        #the temporal dimension
        if verbose: print 'Reducing quantity %s' %quantity
        q_reduced = zeros( number_of_points, Float )

        for k in range(number_of_points):
            q_reduced[k] = reduction( q[:,k] )

        q = q_reduced

    #Post condition: Now q has dimension: number_of_points
    assert len(q.shape) == 1
    assert q.shape[0] == number_of_points


    if verbose:
        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))


    #Create grid and update xll/yll corner and x,y

    #Relative extent
    if easting_min is None:
        xmin = min(x)
    else:
        xmin = easting_min - xllcorner

    if easting_max is None:
        xmax = max(x)
    else:
        xmax = easting_max - xllcorner

    if northing_min is None:
        ymin = min(y)
    else:
        ymin = northing_min - yllcorner

    if northing_max is None:
        ymax = max(y)
    else:
        ymax = northing_max - yllcorner



    if verbose: print 'Creating grid'
    ncols = int((xmax-xmin)/cellsize)+1
    nrows = int((ymax-ymin)/cellsize)+1


    #New absolute reference and coordinates
    newxllcorner = xmin+xllcorner
    newyllcorner = ymin+yllcorner

    x = x+xllcorner-newxllcorner
    y = y+yllcorner-newyllcorner

    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
    assert len(vertex_points.shape) == 2



    grid_points = zeros ( (ncols*nrows, 2), Float )


    for i in xrange(nrows):
        if format.lower() == 'asc':
            yg = i*cellsize
        else:
            #this will flip the order of the y values for ers
            yg = (nrows-i)*cellsize

        for j in xrange(ncols):
            xg = j*cellsize
            k = i*ncols + j

            grid_points[k,0] = xg
            grid_points[k,1] = yg

    #Interpolate
    from least_squares import Interpolation


    #FIXME: This should be done with precrop = True (?), otherwise it'll
    #take forever. With expand_search  set to False, some grid points might
    #miss out.... This will be addressed though Duncan's refactoring of least_squares

    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
                           precrop = False,
                           expand_search = expand_search,
                           verbose = verbose)



    #Interpolate using quantity values
    if verbose: print 'Interpolating'
    grid_values = interp.interpolate(q).flat


    if verbose:
        print 'Interpolated values are in [%f, %f]' %(min(grid_values),
                                                      max(grid_values))

    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
    P = interp.mesh.get_boundary_polygon()
    outside_indices = outside_polygon(grid_points, P, closed=True)

    for i in outside_indices:
        grid_values[i] = NODATA_value




    if format.lower() == 'ers':
        # setup ERS header information
        grid_values = reshape(grid_values,(nrows, ncols))
        header = {}
        header['datum'] = '"' + datum + '"'
        # FIXME The use of hardwired UTM and zone number needs to be made optional
        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
        header['projection'] = '"UTM-' + str(zone) + '"'
        header['coordinatetype'] = 'EN'
        if header['coordinatetype'] == 'LL':
            header['longitude'] = str(newxllcorner)
            header['latitude'] = str(newyllcorner)
        elif header['coordinatetype'] == 'EN':
            header['eastings'] = str(newxllcorner)
            header['northings'] = str(newyllcorner)
        header['nullcellvalue'] = str(NODATA_value)
        header['xdimension'] = str(cellsize)
        header['ydimension'] = str(cellsize)
        header['value'] = '"' + quantity + '"'
        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test


        #Write
        if verbose: print 'Writing %s' %demfile
        import ermapper_grids
        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)

        fid.close()
    else:
        #Write to Ascii format

        #Write prj file
        prjfile = basename_out + '.prj'

        if verbose: print 'Writing %s' %prjfile
        prjid = open(prjfile, 'w')
        prjid.write('Projection    %s\n' %'UTM')
        prjid.write('Zone          %d\n' %zone)
        prjid.write('Datum         %s\n' %datum)
        prjid.write('Zunits        NO\n')
        prjid.write('Units         METERS\n')
        prjid.write('Spheroid      %s\n' %datum)
        prjid.write('Xshift        %d\n' %false_easting)
        prjid.write('Yshift        %d\n' %false_northing)
        prjid.write('Parameters\n')
        prjid.close()



        if verbose: print 'Writing %s' %demfile

        ascid = open(demfile, 'w')

        ascid.write('ncols         %d\n' %ncols)
        ascid.write('nrows         %d\n' %nrows)
        ascid.write('xllcorner     %d\n' %newxllcorner)
        ascid.write('yllcorner     %d\n' %newyllcorner)
        ascid.write('cellsize      %f\n' %cellsize)
        ascid.write('NODATA_value  %d\n' %NODATA_value)


        #Get bounding polygon from mesh
        #P = interp.mesh.get_boundary_polygon()
        #inside_indices = inside_polygon(grid_points, P)

        for i in range(nrows):
            if verbose and i%((nrows+10)/10)==0:
                print 'Doing row %d of %d' %(i, nrows)

            base_index = (nrows-i-1)*ncols

            slice = grid_values[base_index:base_index+ncols]
            s = array2string(slice, max_line_width=sys.maxint)
            ascid.write(s[1:-1] + '\n')


            #print
            #for j in range(ncols):
            #    index = base_index+j#
            #    print grid_values[index],
            #    ascid.write('%f ' %grid_values[index])
            #ascid.write('\n')

        #Close
        ascid.close()
        fid.close()

#Backwards compatibility
def sww2asc(basename_in, basename_out = None,
            quantity = None,
            timestep = None,
            reduction = None,
            cellsize = 10,
            verbose = False,
            origin = None):
    print 'sww2asc will soon be obsoleted - please use sww2dem'
    sww2dem(basename_in,
            basename_out = basename_out,
            quantity = quantity,
            timestep = timestep,
            reduction = reduction,
            cellsize = cellsize,
            verbose = verbose,
            origin = origin,
            datum = 'WGS84',
            format = 'asc')

def sww2ers(basename_in, basename_out = None,
            quantity = None,
            timestep = None,
            reduction = None,
            cellsize = 10,
            verbose = False,
            origin = None,
            datum = 'WGS84'):
    print 'sww2ers will soon be obsoleted - please use sww2dem'
    sww2dem(basename_in,
            basename_out = basename_out,
            quantity = quantity,
            timestep = timestep,
            reduction = reduction,
            cellsize = cellsize,
            verbose = verbose,
            origin = origin,
            datum = datum,
            format = 'ers')
################################# END COMPATIBILITY ##############




def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
                                  verbose=False):
    """Read Digitial Elevation model from the following ASCII format (.asc)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Convert basename_in + '.asc' to NetCDF format (.dem)
    mimicking the ASCII format closely.


    An accompanying file with same basename_in but extension .prj must exist
    and is used to fix the UTM zone, datum, false northings and eastings.

    The prj format is assumed to be as

    Projection    UTM
    Zone          56
    Datum         WGS84
    Zunits        NO
    Units         METERS
    Spheroid      WGS84
    Xshift        0.0000000000
    Yshift        10000000.0000000000
    Parameters
    """

    import os
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, array

    #root, ext = os.path.splitext(basename_in)
    root = basename_in

    ###########################################
    # Read Meta data
    if verbose: print 'Reading METADATA from %s' %root + '.prj'
    metadatafile = open(root + '.prj')
    metalines = metadatafile.readlines()
    metadatafile.close()

    L = metalines[0].strip().split()
    assert L[0].strip().lower() == 'projection'
    projection = L[1].strip()                   #TEXT

    L = metalines[1].strip().split()
    assert L[0].strip().lower() == 'zone'
    zone = int(L[1].strip())

    L = metalines[2].strip().split()
    assert L[0].strip().lower() == 'datum'
    datum = L[1].strip()                        #TEXT

    L = metalines[3].strip().split()
    assert L[0].strip().lower() == 'zunits'     #IGNORE
    zunits = L[1].strip()                       #TEXT

    L = metalines[4].strip().split()
    assert L[0].strip().lower() == 'units'
    units = L[1].strip()                        #TEXT

    L = metalines[5].strip().split()
    assert L[0].strip().lower() == 'spheroid'   #IGNORE
    spheroid = L[1].strip()                     #TEXT

    L = metalines[6].strip().split()
    assert L[0].strip().lower() == 'xshift'
    false_easting = float(L[1].strip())

    L = metalines[7].strip().split()
    assert L[0].strip().lower() == 'yshift'
    false_northing = float(L[1].strip())

    #print false_easting, false_northing, zone, datum


    ###########################################
    #Read DEM data

    datafile = open(basename_in + '.asc')

    if verbose: print 'Reading DEM from %s' %(basename_in + '.asc')
    lines = datafile.readlines()
    datafile.close()

    if verbose: print 'Got', len(lines), ' lines'

    ncols = int(lines[0].split()[1].strip())
    nrows = int(lines[1].split()[1].strip())
    xllcorner = float(lines[2].split()[1].strip())
    yllcorner = float(lines[3].split()[1].strip())
    cellsize = float(lines[4].split()[1].strip())
    NODATA_value = int(lines[5].split()[1].strip())

    assert len(lines) == nrows + 6


    ##########################################


    if basename_out == None:
        netcdfname = root + '.dem'
    else:
        netcdfname = basename_out + '.dem'

    if verbose: print 'Store to NetCDF file %s' %netcdfname
    # NetCDF file definition
    fid = NetCDFFile(netcdfname, 'w')

    #Create new file
    fid.institution = 'Geoscience Australia'
    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
                      'of spatial point data'

    fid.ncols = ncols
    fid.nrows = nrows
    fid.xllcorner = xllcorner
    fid.yllcorner = yllcorner
    fid.cellsize = cellsize
    fid.NODATA_value = NODATA_value

    fid.zone = zone
    fid.false_easting = false_easting
    fid.false_northing = false_northing
    fid.projection = projection
    fid.datum = datum
    fid.units = units


    # dimension definitions
    fid.createDimension('number_of_rows', nrows)
    fid.createDimension('number_of_columns', ncols)

    # variable definitions
    fid.createVariable('elevation', Float, ('number_of_rows',
                                            'number_of_columns'))

    # Get handles to the variables
    elevation = fid.variables['elevation']

    #Store data
    n = len(lines[6:])
    for i, line in enumerate(lines[6:]):
        fields = line.split()
        if verbose and i%((n+10)/10)==0:
            print 'Processing row %d of %d' %(i, nrows)

        elevation[i, :] = array([float(x) for x in fields])

    fid.close()





def ferret2sww(basename_in, basename_out = None,
               verbose = False,
               minlat = None, maxlat = None,
               minlon = None, maxlon = None,
               mint = None, maxt = None, mean_stage = 0,
               origin = None, zscale = 1,
               fail_on_NaN = True,
               NaN_filler = 0,
               elevation = None,
               inverted_bathymetry = False
               ): #FIXME: Bathymetry should be obtained
                                  #from MOST somehow.
                                  #Alternatively from elsewhere
                                  #or, as a last resort,
                                  #specified here.
                                  #The value of -100 will work
                                  #for the Wollongong tsunami
                                  #scenario but is very hacky
    """Convert 'Ferret' NetCDF format for wave propagation to
    sww format native to pyvolution.

    Specify only basename_in and read files of the form
    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
    relative height, x-velocity and y-velocity, respectively.

    Also convert latitude and longitude to UTM. All coordinates are
    assumed to be given in the GDA94 datum.

    min's and max's: If omitted - full extend is used.
    To include a value min may equal it, while max must exceed it.
    Lat and lon are assuemd to be in decimal degrees

    origin is a 3-tuple with geo referenced
    UTM coordinates (zone, easting, northing)

    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
    which means that longitude is the fastest
    varying dimension (row major order, so to speak)

    ferret2sww uses grid points as vertices in a triangular grid
    counting vertices from lower left corner upwards, then right
    """

    import os
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, Int, Int32, searchsorted, zeros, array
    from Numeric import allclose, around

    precision = Float


    #Get NetCDF data
    if verbose: print 'Reading files %s_*.nc' %basename_in
    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)

    if basename_out is None:
        swwname = basename_in + '.sww'
    else:
        swwname = basename_out + '.sww'

    times = file_h.variables['TIME']
    latitudes = file_h.variables['LAT']
    longitudes = file_h.variables['LON']



    #Precision used by most for lat/lon is 4 or 5 decimals
    e_lat = around(file_e.variables['LAT'][:], 5)
    e_lon = around(file_e.variables['LON'][:], 5)

    #Check that files are compatible
    assert allclose(latitudes, file_u.variables['LAT'])
    assert allclose(latitudes, file_v.variables['LAT'])
    assert allclose(latitudes, e_lat)

    assert allclose(longitudes, file_u.variables['LON'])
    assert allclose(longitudes, file_v.variables['LON'])
    assert allclose(longitudes, e_lon)



    if mint == None:
        jmin = 0
    else:
        jmin = searchsorted(times, mint)

    if maxt == None:
        jmax=len(times)
    else:
        jmax = searchsorted(times, maxt)

    if minlat == None:
        kmin=0
    else:
        kmin = searchsorted(latitudes, minlat)

    if maxlat == None:
        kmax = len(latitudes)
    else:
        kmax = searchsorted(latitudes, maxlat)

    if minlon == None:
        lmin=0
    else:
        lmin = searchsorted(longitudes, minlon)

    if maxlon == None:
        lmax = len(longitudes)
    else:
        lmax = searchsorted(longitudes, maxlon)



    times = times[jmin:jmax]
    latitudes = latitudes[kmin:kmax]
    longitudes = longitudes[lmin:lmax]


    if verbose: print 'cropping'
    zname = 'ELEVATION'


    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]

    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
    #        from Numeric import asarray
    #        elevations=elevations.tolist()
    #        elevations.reverse()
    #        elevations=asarray(elevations)
    #    else:
    #        from Numeric import asarray
    #        elevations=elevations.tolist()
    #        elevations.reverse()
    #        elevations=asarray(elevations)
    #        'print hmmm'



    #Get missing values
    nan_ha = file_h.variables['HA'].missing_value[0]
    nan_ua = file_u.variables['UA'].missing_value[0]
    nan_va = file_v.variables['VA'].missing_value[0]
    if hasattr(file_e.variables[zname],'missing_value'):
        nan_e  = file_e.variables[zname].missing_value[0]
    else:
        nan_e = None

    #Cleanup
    from Numeric import sometrue

    missing = (amplitudes == nan_ha)
    if sometrue (missing):
        if fail_on_NaN:
            msg = 'NetCDFFile %s contains missing values'\
                  %(basename_in+'_ha.nc')
            raise msg
        else:
            amplitudes = amplitudes*(missing==0) + missing*NaN_filler

    missing = (uspeed == nan_ua)
    if sometrue (missing):
        if fail_on_NaN:
            msg = 'NetCDFFile %s contains missing values'\
                  %(basename_in+'_ua.nc')
            raise msg
        else:
            uspeed = uspeed*(missing==0) + missing*NaN_filler

    missing = (vspeed == nan_va)
    if sometrue (missing):
        if fail_on_NaN:
            msg = 'NetCDFFile %s contains missing values'\
                  %(basename_in+'_va.nc')
            raise msg
        else:
            vspeed = vspeed*(missing==0) + missing*NaN_filler


    missing = (elevations == nan_e)
    if sometrue (missing):
        if fail_on_NaN:
            msg = 'NetCDFFile %s contains missing values'\
                  %(basename_in+'_e.nc')
            raise msg
        else:
            elevations = elevations*(missing==0) + missing*NaN_filler

    #######



    number_of_times = times.shape[0]
    number_of_latitudes = latitudes.shape[0]
    number_of_longitudes = longitudes.shape[0]

    assert amplitudes.shape[0] == number_of_times
    assert amplitudes.shape[1] == number_of_latitudes
    assert amplitudes.shape[2] == number_of_longitudes

    if verbose:
        print '------------------------------------------------'
        print 'Statistics:'
        print '  Extent (lat/lon):'
        print '    lat in [%f, %f], len(lat) == %d'\
              %(min(latitudes.flat), max(latitudes.flat),
                len(latitudes.flat))
        print '    lon in [%f, %f], len(lon) == %d'\
              %(min(longitudes.flat), max(longitudes.flat),
                len(longitudes.flat))
        print '    t in [%f, %f], len(t) == %d'\
              %(min(times.flat), max(times.flat), len(times.flat))

        q = amplitudes.flat
        name = 'Amplitudes (ha) [cm]'
        print '  %s in [%f, %f]' %(name, min(q), max(q))

        q = uspeed.flat
        name = 'Speeds (ua) [cm/s]'
        print '  %s in [%f, %f]' %(name, min(q), max(q))

        q = vspeed.flat
        name = 'Speeds (va) [cm/s]'
        print '  %s in [%f, %f]' %(name, min(q), max(q))

        q = elevations.flat
        name = 'Elevations (e) [m]'
        print '  %s in [%f, %f]' %(name, min(q), max(q))


    #print number_of_latitudes, number_of_longitudes
    number_of_points = number_of_latitudes*number_of_longitudes
    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2


    file_h.close()
    file_u.close()
    file_v.close()
    file_e.close()


    # NetCDF file definition
    outfile = NetCDFFile(swwname, 'w')

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\
                          %(basename_in + '_ha.nc',
                            basename_in + '_ua.nc',
                            basename_in + '_va.nc',
                            basename_in + '_e.nc')


    #For sww compatibility
    outfile.smoothing = 'Yes'
    outfile.order = 1

    #Start time in seconds since the epoch (midnight 1/1/1970)
    outfile.starttime = starttime = times[0]
    times = times - starttime  #Store relative times

    # dimension definitions
    outfile.createDimension('number_of_volumes', number_of_volumes)

    outfile.createDimension('number_of_vertices', 3)
    outfile.createDimension('number_of_points', number_of_points)


    #outfile.createDimension('number_of_timesteps', len(times))
    outfile.createDimension('number_of_timesteps', len(times))

    # variable definitions
    outfile.createVariable('x', precision, ('number_of_points',))
    outfile.createVariable('y', precision, ('number_of_points',))
    outfile.createVariable('elevation', precision, ('number_of_points',))

    #FIXME: Backwards compatibility
    outfile.createVariable('z', precision, ('number_of_points',))
    #################################

    outfile.createVariable('volumes', Int, ('number_of_volumes',
                                            'number_of_vertices'))

    outfile.createVariable('time', precision,
                           ('number_of_timesteps',))

    outfile.createVariable('stage', precision,
                           ('number_of_timesteps',
                            'number_of_points'))

    outfile.createVariable('xmomentum', precision,
                           ('number_of_timesteps',
                            'number_of_points'))

    outfile.createVariable('ymomentum', precision,
                           ('number_of_timesteps',
                            'number_of_points'))


    #Store
    from coordinate_transforms.redfearn import redfearn
    x = zeros(number_of_points, Float)  #Easting
    y = zeros(number_of_points, Float)  #Northing


    if verbose: print 'Making triangular grid'
    #Check zone boundaries
    refzone, _, _ = redfearn(latitudes[0],longitudes[0])

    vertices = {}
    i = 0
    for k, lat in enumerate(latitudes):       #Y direction
        for l, lon in enumerate(longitudes):  #X direction

            vertices[l,k] = i

            zone, easting, northing = redfearn(lat,lon)

            msg = 'Zone boundary crossed at longitude =', lon
            #assert zone == refzone, msg
            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
            x[i] = easting
            y[i] = northing
            i += 1


    #Construct 2 triangles per 'rectangular' element
    volumes = []
    for l in range(number_of_longitudes-1):    #X direction
        for k in range(number_of_latitudes-1): #Y direction
            v1 = vertices[l,k+1]
            v2 = vertices[l,k]
            v3 = vertices[l+1,k+1]
            v4 = vertices[l+1,k]

            volumes.append([v1,v2,v3]) #Upper element
            volumes.append([v4,v3,v2]) #Lower element

    volumes = array(volumes)

    if origin == None:
        zone = refzone
        xllcorner = min(x)
        yllcorner = min(y)
    else:
        zone = origin[0]
        xllcorner = origin[1]
        yllcorner = origin[2]


    outfile.xllcorner = xllcorner
    outfile.yllcorner = yllcorner
    outfile.zone = zone


    if elevation is not None:
        z = elevation
    else:
        if inverted_bathymetry:
            z = -1*elevations
        else:
            z = elevations
        #FIXME: z should be obtained from MOST and passed in here

    from Numeric import resize
    z = resize(z,outfile.variables['z'][:].shape)
    outfile.variables['x'][:] = x - xllcorner
    outfile.variables['y'][:] = y - yllcorner
    outfile.variables['z'][:] = z             #FIXME HACK
    outfile.variables['elevation'][:] = z
    outfile.variables['time'][:] = times   #Store time relative
    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64



    #Time stepping
    stage = outfile.variables['stage']
    xmomentum = outfile.variables['xmomentum']
    ymomentum = outfile.variables['ymomentum']

    if verbose: print 'Converting quantities'
    n = len(times)
    for j in range(n):
        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
        i = 0
        for k in range(number_of_latitudes):      #Y direction
            for l in range(number_of_longitudes): #X direction
                w = zscale*amplitudes[j,k,l]/100 + mean_stage
                stage[j,i] = w
                h = w - z[i]
                xmomentum[j,i] = uspeed[j,k,l]/100*h
                ymomentum[j,i] = vspeed[j,k,l]/100*h
                i += 1

    #outfile.close()

    #FIXME: Refactor using code from file_function.statistics
    #Something like print swwstats(swwname)
    if verbose:
        x = outfile.variables['x'][:]
        y = outfile.variables['y'][:]
        times = outfile.variables['time'][:]
        print '------------------------------------------------'
        print 'Statistics of output file:'
        print '  Name: %s' %swwname
        print '  Reference:'
        print '    Lower left corner: [%f, %f]'\
              %(xllcorner, yllcorner)
        print '    Start time: %f' %starttime
        print '  Extent:'
        print '    x [m] in [%f, %f], len(x) == %d'\
              %(min(x.flat), max(x.flat), len(x.flat))
        print '    y [m] in [%f, %f], len(y) == %d'\
              %(min(y.flat), max(y.flat), len(y.flat))
        print '    t [s] in [%f, %f], len(t) == %d'\
              %(min(times), max(times), len(times))
        print '  Quantities [SI units]:'
        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
            q = outfile.variables[name][:].flat
            print '    %s in [%f, %f]' %(name, min(q), max(q))



    outfile.close()





def timefile2netcdf(filename, quantity_names = None):
    """Template for converting typical text files with time series to
    NetCDF tms file.


    The file format is assumed to be either two fields separated by a comma:

        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...

    E.g

      31/08/04 00:00:00, 1.328223 0 0
      31/08/04 00:15:00, 1.292912 0 0

    will provide a time dependent function f(t) with three attributes

    filename is assumed to be the rootname with extenisons .txt and .sww
    """

    import time, calendar
    from Numeric import array
    from config import time_format
    from util import ensure_numeric


    fid = open(filename + '.txt')
    line = fid.readline()
    fid.close()

    fields = line.split(',')
    msg = 'File %s must have the format date, value0 value1 value2 ...'
    assert len(fields) == 2, msg

    try:
        starttime = calendar.timegm(time.strptime(fields[0], time_format))
    except ValueError:
        msg = 'First field in file %s must be' %filename
        msg += ' date-time with format %s.\n' %time_format
        msg += 'I got %s instead.' %fields[0]
        raise msg


    #Split values
    values = []
    for value in fields[1].split():
        values.append(float(value))

    q = ensure_numeric(values)

    msg = 'ERROR: File must contain at least one independent value'
    assert len(q.shape) == 1, msg



    #Read times proper
    from Numeric import zeros, Float, alltrue
    from config import time_format
    import time, calendar

    fid = open(filename + '.txt')
    lines = fid.readlines()
    fid.close()

    N = len(lines)
    d = len(q)

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

    for i, line in enumerate(lines):
        fields = line.split(',')
        realtime = calendar.timegm(time.strptime(fields[0], time_format))

        T[i] = realtime - starttime

        for j, value in enumerate(fields[1].split()):
            Q[i, j] = float(value)

    msg = 'File %s must list time as a monotonuosly ' %filename
    msg += 'increasing sequence'
    assert alltrue( T[1:] - T[:-1] > 0 ), msg


    #Create NetCDF file
    from Scientific.IO.NetCDF import NetCDFFile

    fid = NetCDFFile(filename + '.tms', 'w')


    fid.institution = 'Geoscience Australia'
    fid.description = 'Time series'


    #Reference point
    #Start time in seconds since the epoch (midnight 1/1/1970)
    #FIXME: Use Georef
    fid.starttime = starttime


    # dimension definitions
    #fid.createDimension('number_of_volumes', self.number_of_volumes)
    #fid.createDimension('number_of_vertices', 3)


    fid.createDimension('number_of_timesteps', len(T))

    fid.createVariable('time', Float, ('number_of_timesteps',))

    fid.variables['time'][:] = T

    for i in range(Q.shape[1]):
        try:
            name = quantity_names[i]
        except:
            name = 'Attribute%d'%i

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

    fid.close()


def extent_sww(file_name):
    """
    Read in an sww file.

    Input;
    file_name - the sww file

    Output;
    z - Vector of bed elevation
    volumes - Array.  Each row has 3 values, representing
    the vertices that define the volume
    time - Vector of the times where there is stage information
    stage - array with respect to time and vertices (x,y)
    """


    from Scientific.IO.NetCDF import NetCDFFile

    #Check contents
    #Get NetCDF
    fid = NetCDFFile(file_name, 'r')

    # Get the variables
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    stage = fid.variables['stage'][:]
    #print "stage",stage
    #print "stage.shap",stage.shape
    #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat)
    #print "min(stage)",min(stage)

    fid.close()

    return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]


def sww2domain(filename,boundary=None,t=None,\
               fail_if_NaN=True,NaN_filler=0\
               ,verbose = False,very_verbose = False):
    """
    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))

    Boundary is not recommended if domain.smooth is not selected, as it
    uses unique coordinates, but not unique boundaries. This means that
    the boundary file will not be compatable with the coordinates, and will
    give a different final boundary, or crash.
    """
    NaN=9.969209968386869e+036
    #initialise NaN.

    from Scientific.IO.NetCDF import NetCDFFile
    from shallow_water import Domain
    from Numeric import asarray, transpose, resize

    if verbose: print 'Reading from ', filename
    fid = NetCDFFile(filename, 'r')    #Open existing file for read
    time = fid.variables['time']       #Timesteps
    if t is None:
        t = time[-1]
    time_interp = get_time_interp(time,t)

    # Get the variables as Numeric arrays
    x = fid.variables['x'][:]             #x-coordinates of vertices
    y = fid.variables['y'][:]             #y-coordinates of vertices
    elevation = fid.variables['elevation']     #Elevation
    stage = fid.variables['stage']     #Water level
    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction

    starttime = fid.starttime[0]
    volumes = fid.variables['volumes'][:] #Connectivity
    coordinates=transpose(asarray([x.tolist(),y.tolist()]))

    conserved_quantities = []
    interpolated_quantities = {}
    other_quantities = []

    # get geo_reference
    #sww files don't have to have a geo_ref
    try:
        geo_reference = Geo_reference(NetCDFObject=fid)
    except: #AttributeError, e:
        geo_reference = None

    if verbose: print '    getting quantities'
    for quantity in fid.variables.keys():
        dimensions = fid.variables[quantity].dimensions
        if 'number_of_timesteps' in dimensions:
            conserved_quantities.append(quantity)
            interpolated_quantities[quantity]=\
                  interpolated_quantity(fid.variables[quantity][:],time_interp)
        else: other_quantities.append(quantity)

    other_quantities.remove('x')
    other_quantities.remove('y')
    other_quantities.remove('z')
    other_quantities.remove('volumes')

    conserved_quantities.remove('time')

    if verbose: print '    building domain'
    #    From domain.Domain:
    #    domain = Domain(coordinates, volumes,\
    #                    conserved_quantities = conserved_quantities,\
    #                    other_quantities = other_quantities,zone=zone,\
    #                    xllcorner=xllcorner, yllcorner=yllcorner)

    #   From shallow_water.Domain:
    coordinates=coordinates.tolist()
    volumes=volumes.tolist()
    #FIXME:should this be in mesh?(peter row)
    if fid.smoothing == 'Yes': unique = False
    else: unique = True
    if unique:
        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)


    try:    
        domain = Domain(coordinates, volumes, boundary)
    except AssertionError, e:
        fid.close()
        msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
        raise msg

    if not boundary is None:
        domain.boundary = boundary

    domain.geo_reference = geo_reference

    domain.starttime=float(starttime)+float(t)
    domain.time=0.0

    for quantity in other_quantities:
        try:
            NaN = fid.variables[quantity].missing_value
        except:
            pass #quantity has no missing_value number
        X = fid.variables[quantity][:]
        if very_verbose:
            print '       ',quantity
            print '        NaN =',NaN
            print '        max(X)'
            print '       ',max(X)
            print '        max(X)==NaN'
            print '       ',max(X)==NaN
            print ''
        if (max(X)==NaN) or (min(X)==NaN):
            if fail_if_NaN:
                msg = 'quantity "%s" contains no_data entry'%quantity
                raise msg
            else:
                data = (X<>NaN)
                X = (X*data)+(data==0)*NaN_filler
        if unique:
            X = resize(X,(len(X)/3,3))
        domain.set_quantity(quantity,X)
    #
    for quantity in conserved_quantities:
        try:
            NaN = fid.variables[quantity].missing_value
        except:
            pass #quantity has no missing_value number
        X = interpolated_quantities[quantity]
        if very_verbose:
            print '       ',quantity
            print '        NaN =',NaN
            print '        max(X)'
            print '       ',max(X)
            print '        max(X)==NaN'
            print '       ',max(X)==NaN
            print ''
        if (max(X)==NaN) or (min(X)==NaN):
            if fail_if_NaN:
                msg = 'quantity "%s" contains no_data entry'%quantity
                raise msg
            else:
                data = (X<>NaN)
                X = (X*data)+(data==0)*NaN_filler
        if unique:
            X = resize(X,(X.shape[0]/3,3))
        domain.set_quantity(quantity,X)
        
    fid.close()
    return domain

def interpolated_quantity(saved_quantity,time_interp):

    #given an index and ratio, interpolate quantity with respect to time.
    index,ratio = time_interp
    Q = saved_quantity
    if ratio > 0:
        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
    else:
        q = Q[index]
    #Return vector of interpolated values
    return q

def get_time_interp(time,t=None):
    #Finds the ratio and index for time interpolation.
    #It is borrowed from previous pyvolution code.
    if t is None:
        t=time[-1]
        index = -1
        ratio = 0.
    else:
        T = time
        tau = t
        index=0
        msg = 'Time interval derived from file %s [%s:%s]'\
            %('FIXMEfilename', T[0], T[-1])
        msg += ' does not match model time: %s' %tau
        if tau < time[0]: raise msg
        if tau > time[-1]: raise msg
        while tau > time[index]: index += 1
        while tau < time[index]: index -= 1
        if tau == time[index]:
            #Protect against case where tau == time[-1] (last time)
            # - also works in general when tau == time[i]
            ratio = 0
        else:
            #t is now between index and index+1
            ratio = (tau - time[index])/(time[index+1] - time[index])
    return (index,ratio)


def weed(coordinates,volumes,boundary = None):
    if type(coordinates)=='array':
        coordinates = coordinates.tolist()
    if type(volumes)=='array':
        volumes = volumes.tolist()

    unique = False
    point_dict = {}
    same_point = {}
    for i in range(len(coordinates)):
        point = tuple(coordinates[i])
        if point_dict.has_key(point):
            unique = True
            same_point[i]=point
            #to change all point i references to point j
        else:
            point_dict[point]=i
            same_point[i]=point

    coordinates = []
    i = 0
    for point in point_dict.keys():
        point = tuple(point)
        coordinates.append(list(point))
        point_dict[point]=i
        i+=1


    for volume in volumes:
        for i in range(len(volume)):
            index = volume[i]
            if index>-1:
                volume[i]=point_dict[same_point[index]]

    new_boundary = {}
    if not boundary is None:
        for segment in boundary.keys():
            point0 = point_dict[same_point[segment[0]]]
            point1 = point_dict[same_point[segment[1]]]
            label = boundary[segment]
            #FIXME should the bounday attributes be concaterated
            #('exterior, pond') or replaced ('pond')(peter row)

            if new_boundary.has_key((point0,point1)):
                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
                                              #+','+label

            elif new_boundary.has_key((point1,point0)):
                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
                                              #+','+label
            else: new_boundary[(point0,point1)]=label

        boundary = new_boundary

    return coordinates,volumes,boundary










###############################################
#OBSOLETE STUFF
#Native checkpoint format.
#Information needed to recreate a state is preserved
#FIXME: Rethink and maybe use netcdf format
def cpt_variable_writer(filename, t, v0, v1, v2):
    """Store all conserved quantities to file
    """

    M, N = v0.shape

    FN = create_filename(filename, 'cpt', M, t)
    #print 'Writing to %s' %FN

    fid = open(FN, 'w')
    for i in range(M):
        for j in range(N):
            fid.write('%.16e ' %v0[i,j])
        for j in range(N):
            fid.write('%.16e ' %v1[i,j])
        for j in range(N):
            fid.write('%.16e ' %v2[i,j])

        fid.write('\n')
    fid.close()


def cpt_variable_reader(filename, t, v0, v1, v2):
    """Store all conserved quantities to file
    """

    M, N = v0.shape

    FN = create_filename(filename, 'cpt', M, t)
    #print 'Reading from %s' %FN

    fid = open(FN)


    for i in range(M):
        values = fid.readline().split() #Get one line

        for j in range(N):
            v0[i,j] = float(values[j])
            v1[i,j] = float(values[3+j])
            v2[i,j] = float(values[6+j])

    fid.close()

def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
    """Writes x,y,z,z,z coordinates of triangles constituting the bed
    elevation.
    FIXME: Not in use pt
    """

    M, N = v0.shape


    print X0
    import sys; sys.exit()
    FN = create_filename(filename, 'cpt', M)
    print 'Writing to %s' %FN

    fid = open(FN, 'w')
    for i in range(M):
        for j in range(2):
            fid.write('%.16e ' %X0[i,j])   #x, y
        for j in range(N):
            fid.write('%.16e ' %v0[i,j])       #z,z,z,

        for j in range(2):
            fid.write('%.16e ' %X1[i,j])   #x, y
        for j in range(N):
            fid.write('%.16e ' %v1[i,j])

        for j in range(2):
            fid.write('%.16e ' %X2[i,j])   #x, y
        for j in range(N):
            fid.write('%.16e ' %v2[i,j])

        fid.write('\n')
    fid.close()



#Function for storing out to e.g. visualisation
#FIXME: Do we want this?
#FIXME: Not done yet for this version
def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
    """Writes x,y,z coordinates of triangles constituting the bed elevation.
    """

    M, N = v0.shape

    FN = create_filename(filename, 'dat', M)
    #print 'Writing to %s' %FN

    fid = open(FN, 'w')
    for i in range(M):
        for j in range(2):
            fid.write('%f ' %X0[i,j])   #x, y
        fid.write('%f ' %v0[i,0])       #z

        for j in range(2):
            fid.write('%f ' %X1[i,j])   #x, y
        fid.write('%f ' %v1[i,0])       #z

        for j in range(2):
            fid.write('%f ' %X2[i,j])   #x, y
        fid.write('%f ' %v2[i,0])       #z

        fid.write('\n')
    fid.close()



def dat_variable_writer(filename, t, v0, v1, v2):
    """Store water height to file
    """

    M, N = v0.shape

    FN = create_filename(filename, 'dat', M, t)
    #print 'Writing to %s' %FN

    fid = open(FN, 'w')
    for i in range(M):
        fid.write('%.4f ' %v0[i,0])
        fid.write('%.4f ' %v1[i,0])
        fid.write('%.4f ' %v2[i,0])

        fid.write('\n')
    fid.close()


def read_sww(filename):
    """Read sww Net CDF file containing Shallow Water Wave simulation

    The integer array volumes is of shape Nx3 where N is the number of
    triangles in the mesh.

    Each entry in volumes is an index into the x,y arrays (the location).

    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
    number_of_timesteps, number_of_points.

    The momentum is not always stored.

    """
    from Scientific.IO.NetCDF import NetCDFFile
    print 'Reading from ', filename
    fid = NetCDFFile(filename, 'r')    #Open existing file for read
#latitude, longitude
    # Get the variables as Numeric arrays
    x = fid.variables['x']             #x-coordinates of vertices
    y = fid.variables['y']             #y-coordinates of vertices
    z = fid.variables['elevation']     #Elevation
    time = fid.variables['time']       #Timesteps
    stage = fid.variables['stage']     #Water level
    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction

    volumes = fid.variables['volumes'] #Connectivity

    #FIXME (Ole): What is this?
    #             Why isn't anything returned?
    #             Where's the unit test?

def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
                 verbose=False):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
    """

    import os
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import Float, zeros, sum, reshape, equal

    root = basename_in
    inname = root + '.dem'

    #Open existing netcdf file to read
    infile = NetCDFFile(inname, 'r')
    if verbose: print 'Reading DEM from %s' %inname

    #Read metadata
    ncols = infile.ncols[0]
    nrows = infile.nrows[0]
    xllcorner = infile.xllcorner[0]
    yllcorner = infile.yllcorner[0]
    cellsize = infile.cellsize[0]
    NODATA_value = infile.NODATA_value[0]
    zone = infile.zone[0]
    false_easting = infile.false_easting[0]
    false_northing = infile.false_northing[0]
    projection = infile.projection
    datum = infile.datum
    units = infile.units

    dem_elevation = infile.variables['elevation']

    #Get output file name
    if basename_out == None:
        outname = root + '_' + repr(cellsize_new) + '.dem'
    else:
        outname = basename_out + '.dem'

    if verbose: print 'Write decimated NetCDF file to %s' %outname

    #Determine some dimensions for decimated grid
    (nrows_stencil, ncols_stencil) = stencil.shape
    x_offset = ncols_stencil / 2
    y_offset = nrows_stencil / 2
    cellsize_ratio = int(cellsize_new / cellsize)
    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio

    #Open netcdf file for output
    outfile = NetCDFFile(outname, 'w')

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
                           'of spatial point data'
    #Georeferencing
    outfile.zone = zone
    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units

    outfile.cellsize = cellsize_new
    outfile.NODATA_value = NODATA_value
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.xllcorner = xllcorner + (x_offset * cellsize)
    outfile.yllcorner = yllcorner + (y_offset * cellsize)
    outfile.ncols = ncols_new
    outfile.nrows = nrows_new

    # dimension definition
    outfile.createDimension('number_of_points', nrows_new*ncols_new)

    # variable definition
    outfile.createVariable('elevation', Float, ('number_of_points',))

    # Get handle to the variable
    elevation = outfile.variables['elevation']

    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))

    #Store data
    global_index = 0
    for i in range(nrows_new):
        if verbose: print 'Processing row %d of %d' %(i, nrows_new)
        lower_index = global_index
        telev =  zeros(ncols_new, Float)
        local_index = 0
        trow = i * cellsize_ratio

        for j in range(ncols_new):
            tcol = j * cellsize_ratio
            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]

            #if dem contains 1 or more NODATA_values set value in
            #decimated dem to NODATA_value, else compute decimated
            #value using stencil
            if sum(sum(equal(tmp, NODATA_value))) > 0:
                telev[local_index] = NODATA_value
            else:
                telev[local_index] = sum(sum(tmp * stencil))

            global_index += 1
            local_index += 1

        upper_index = global_index

        elevation[lower_index:upper_index] = telev

    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'

    infile.close()
    outfile.close()



def sww2asc_obsolete(basename_in, basename_out = None,
            quantity = None,
            timestep = None,
            reduction = None,
            cellsize = 10,
            verbose = False,
            origin = None):
    """Read SWW file and convert to Digitial Elevation model format (.asc)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Also write accompanying file with same basename_in but extension .prj
    used to fix the UTM zone, datum, false northings and eastings.

    The prj format is assumed to be as

    Projection    UTM
    Zone          56
    Datum         WGS84
    Zunits        NO
    Units         METERS
    Spheroid      WGS84
    Xshift        0.0000000000
    Yshift        10000000.0000000000
    Parameters


    if quantity is given, out values from quantity otherwise default to
    elevation

    if timestep (an index) is given, output quantity at that timestep

    if reduction is given use that to reduce quantity over all timesteps.

    """
    from Numeric import array, Float, concatenate, NewAxis, zeros,\
         sometrue
    from utilities.polygon import inside_polygon

    #FIXME: Should be variable
    datum = 'WGS84'
    false_easting = 500000
    false_northing = 10000000

    if quantity is None:
        quantity = 'elevation'

    if reduction is None:
        reduction = max

    if basename_out is None:
        basename_out = basename_in + '_%s' %quantity

    swwfile = basename_in + '.sww'
    ascfile = basename_out + '.asc'
    prjfile = basename_out + '.prj'


    if verbose: print 'Reading from %s' %swwfile
    #Read sww file
    from Scientific.IO.NetCDF import NetCDFFile
    fid = NetCDFFile(swwfile)

    #Get extent and reference
    x = fid.variables['x'][:]
    y = fid.variables['y'][:]
    volumes = fid.variables['volumes'][:]

    ymin = min(y); ymax = max(y)
    xmin = min(x); xmax = max(x)

    number_of_timesteps = fid.dimensions['number_of_timesteps']
    number_of_points = fid.dimensions['number_of_points']
    if origin is None:

        #Get geo_reference
        #sww files don't have to have a geo_ref
        try:
            geo_reference = Geo_reference(NetCDFObject=fid)
        except AttributeError, e:
            geo_reference = Geo_reference() #Default georef object

        xllcorner = geo_reference.get_xllcorner()
        yllcorner = geo_reference.get_yllcorner()
        zone = geo_reference.get_zone()
    else:
        zone = origin[0]
        xllcorner = origin[1]
        yllcorner = origin[2]


    #Get quantity and reduce if applicable
    if verbose: print 'Reading quantity %s' %quantity

    if quantity.lower() == 'depth':
        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
    else:
        q = fid.variables[quantity][:]


    if len(q.shape) == 2:
        if verbose: print 'Reducing quantity %s' %quantity
        q_reduced = zeros( number_of_points, Float )

        for k in range(number_of_points):
            q_reduced[k] = reduction( q[:,k] )

        q = q_reduced

    #Now q has dimension: number_of_points

    #Create grid and update xll/yll corner and x,y
    if verbose: print 'Creating grid'
    ncols = int((xmax-xmin)/cellsize)+1
    nrows = int((ymax-ymin)/cellsize)+1

    newxllcorner = xmin+xllcorner
    newyllcorner = ymin+yllcorner

    x = x+xllcorner-newxllcorner
    y = y+yllcorner-newyllcorner

    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
    assert len(vertex_points.shape) == 2


    from Numeric import zeros, Float
    grid_points = zeros ( (ncols*nrows, 2), Float )


    for i in xrange(nrows):
        yg = i*cellsize
        for j in xrange(ncols):
            xg = j*cellsize
            k = i*ncols + j

            grid_points[k,0] = xg
            grid_points[k,1] = yg

    #Interpolate
    from least_squares import Interpolation


    #FIXME: This should be done with precrop = True, otherwise it'll
    #take forever. With expand_search  set to False, some grid points might
    #miss out....

    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
                           precrop = False, expand_search = True,
                           verbose = verbose)

    #Interpolate using quantity values
    if verbose: print 'Interpolating'
    grid_values = interp.interpolate(q).flat

    #Write
    #Write prj file
    if verbose: print 'Writing %s' %prjfile
    prjid = open(prjfile, 'w')
    prjid.write('Projection    %s\n' %'UTM')
    prjid.write('Zone          %d\n' %zone)
    prjid.write('Datum         %s\n' %datum)
    prjid.write('Zunits        NO\n')
    prjid.write('Units         METERS\n')
    prjid.write('Spheroid      %s\n' %datum)
    prjid.write('Xshift        %d\n' %false_easting)
    prjid.write('Yshift        %d\n' %false_northing)
    prjid.write('Parameters\n')
    prjid.close()



    if verbose: print 'Writing %s' %ascfile
    NODATA_value = -9999

    ascid = open(ascfile, 'w')

    ascid.write('ncols         %d\n' %ncols)
    ascid.write('nrows         %d\n' %nrows)
    ascid.write('xllcorner     %d\n' %newxllcorner)
    ascid.write('yllcorner     %d\n' %newyllcorner)
    ascid.write('cellsize      %f\n' %cellsize)
    ascid.write('NODATA_value  %d\n' %NODATA_value)


    #Get bounding polygon from mesh
    P = interp.mesh.get_boundary_polygon()
    inside_indices = inside_polygon(grid_points, P)

    for i in range(nrows):
        if verbose and i%((nrows+10)/10)==0:
            print 'Doing row %d of %d' %(i, nrows)

        for j in range(ncols):
            index = (nrows-i-1)*ncols+j

            if sometrue(inside_indices == index):
                ascid.write('%f ' %grid_values[index])
            else:
                ascid.write('%d ' %NODATA_value)

        ascid.write('\n')

    #Close
    ascid.close()
    fid.close()

#********************
#*** END OF OBSOLETE FUNCTIONS
#***************


def tsh2sww(filename, verbose=False): #test_tsh2sww
    """
    to check if a tsh/msh file 'looks' good.
    """

    #FIXME Move to data_manager
    from shallow_water import Domain
    from pmesh2domain import pmesh_to_domain_instance
    import time, os
    from data_manager import get_dataobject
    from os import sep, path
    from util import mean

    if verbose == True:print 'Creating domain from', filename
    domain = pmesh_to_domain_instance(filename, Domain)
    if verbose == True:print "Number of triangles = ", len(domain)

    domain.smooth = True
    domain.format = 'sww'   #Native netcdf visualisation format
    file_path, filename = path.split(filename)
    filename, ext = path.splitext(filename)
    domain.filename = filename
    domain.reduction = mean
    if verbose == True:print "file_path",file_path
    if file_path == "":file_path = "."
    domain.set_datadir(file_path)

    if verbose == True:
        print "Output written to " + domain.get_datadir() + sep + \
              domain.filename + "." + domain.format
    sww = get_dataobject(domain)
    sww.store_connectivity()
    sww.store_timestep('stage')


def asc_csiro2sww(bath_dir,
                  elevation_dir,
                  ucur_dir,
                  vcur_dir,
                  sww_file,
                  minlat = None, maxlat = None,
                  minlon = None, maxlon = None,
                  zscale=1,
                  mean_stage = 0,
                  fail_on_NaN = True,
                  elevation_NaN_filler = 0,
                  bath_prefix='ba',
                  elevation_prefix='el',
                  verbose=False):
    """
    Produce an sww boundary file, from esri ascii data from CSIRO.

    Also convert latitude and longitude to UTM. All coordinates are
    assumed to be given in the GDA94 datum.

    assume:
    All files are in esri ascii format

    4 types of information
    bathymetry
    elevation
    u velocity
    v velocity

    Assumptions
    The metadata of all the files is the same
    Each type is in a seperate directory
    One bath file with extention .000
    The time period is less than 24hrs and uniform.
    """
    from Scientific.IO.NetCDF import NetCDFFile

    from coordinate_transforms.redfearn import redfearn

    precision = Float # So if we want to change the precision its done here

    # go in to the bath dir and load the only file,
    bath_files = os.listdir(bath_dir)
    #print "bath_files",bath_files

    #fixme: make more general?
    bath_file = bath_files[0]
    bath_dir_file =  bath_dir + os.sep + bath_file
    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
    #print "bath_metadata",bath_metadata
    #print "bath_grid",bath_grid

    #Use the date.time of the bath file as a basis for
    #the start time for other files
    base_start = bath_file[-12:]

    #go into the elevation dir and load the 000 file
    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
                         + base_start
    #elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
    #print "elevation_dir_file",elevation_dir_file
    #print "elevation_grid", elevation_grid

    elevation_files = os.listdir(elevation_dir)
    ucur_files = os.listdir(ucur_dir)
    vcur_files = os.listdir(vcur_dir)

    # the first elevation file should be the
    # file with the same base name as the bath data
    #print "elevation_files[0]",elevation_files[0]
    #print "'el' + base_start",'el' + base_start
    assert elevation_files[0] == 'el' + base_start

    #assert bath_metadata == elevation_metadata



    number_of_latitudes = bath_grid.shape[0]
    number_of_longitudes = bath_grid.shape[1]
    #number_of_times = len(os.listdir(elevation_dir))
    #number_of_points = number_of_latitudes*number_of_longitudes
    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2

    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
                  for x in range(number_of_longitudes)]
    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
                 for y in range(number_of_latitudes)]

     # reverse order of lat, so the fist lat represents the first grid row
    latitudes.reverse()

    #print "latitudes - before _get_min_max_indexes",latitudes
    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes,longitudes,
                                                 minlat=minlat, maxlat=maxlat,
                                                 minlon=minlon, maxlon=maxlon)


    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
    #print "bath_grid",bath_grid
    latitudes = latitudes[kmin:kmax]
    longitudes = longitudes[lmin:lmax]
    number_of_latitudes = len(latitudes)
    number_of_longitudes = len(longitudes)
    number_of_times = len(os.listdir(elevation_dir))
    number_of_points = number_of_latitudes*number_of_longitudes
    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
    #print "number_of_points",number_of_points

    #Work out the times
    if len(elevation_files) > 1:
        # Assume: The time period is less than 24hrs.
        time_period = (int(elevation_files[1][-3:]) - \
                      int(elevation_files[0][-3:]))*60*60
        times = [x*time_period for x in range(len(elevation_files))]
    else:
        times = [0.0]
    #print "times", times
    #print "number_of_latitudes", number_of_latitudes
    #print "number_of_longitudes", number_of_longitudes
    #print "number_of_times", number_of_times
    #print "latitudes", latitudes
    #print "longitudes", longitudes


    if verbose:
        print '------------------------------------------------'
        print 'Statistics:'
        print '  Extent (lat/lon):'
        print '    lat in [%f, %f], len(lat) == %d'\
              %(min(latitudes), max(latitudes),
                len(latitudes))
        print '    lon in [%f, %f], len(lon) == %d'\
              %(min(longitudes), max(longitudes),
                len(longitudes))
        print '    t in [%f, %f], len(t) == %d'\
              %(min(times), max(times), len(times))

    ######### WRITE THE SWW FILE #############
    # NetCDF file definition
    outfile = NetCDFFile(sww_file, 'w')

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'Converted from XXX'


    #For sww compatibility
    outfile.smoothing = 'Yes'
    outfile.order = 1

    #Start time in seconds since the epoch (midnight 1/1/1970)
    outfile.starttime = starttime = times[0]


    # dimension definitions
    outfile.createDimension('number_of_volumes', number_of_volumes)

    outfile.createDimension('number_of_vertices', 3)
    outfile.createDimension('number_of_points', number_of_points)
    outfile.createDimension('number_of_timesteps', number_of_times)

    # variable definitions
    outfile.createVariable('x', precision, ('number_of_points',))
    outfile.createVariable('y', precision, ('number_of_points',))
    outfile.createVariable('elevation', precision, ('number_of_points',))

    #FIXME: Backwards compatibility
    outfile.createVariable('z', precision, ('number_of_points',))
    #################################

    outfile.createVariable('volumes', Int, ('number_of_volumes',
                                            'number_of_vertices'))

    outfile.createVariable('time', precision,
                           ('number_of_timesteps',))

    outfile.createVariable('stage', precision,
                           ('number_of_timesteps',
                            'number_of_points'))

    outfile.createVariable('xmomentum', precision,
                           ('number_of_timesteps',
                            'number_of_points'))

    outfile.createVariable('ymomentum', precision,
                           ('number_of_timesteps',
                            'number_of_points'))

    #Store
    from coordinate_transforms.redfearn import redfearn
    x = zeros(number_of_points, Float)  #Easting
    y = zeros(number_of_points, Float)  #Northing

    if verbose: print 'Making triangular grid'
    #Get zone of 1st point.
    refzone, _, _ = redfearn(latitudes[0],longitudes[0])

    vertices = {}
    i = 0
    for k, lat in enumerate(latitudes):
        for l, lon in enumerate(longitudes):

            vertices[l,k] = i

            zone, easting, northing = redfearn(lat,lon)

            msg = 'Zone boundary crossed at longitude =', lon
            #assert zone == refzone, msg
            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
            x[i] = easting
            y[i] = northing
            i += 1


    #Construct 2 triangles per 'rectangular' element
    volumes = []
    for l in range(number_of_longitudes-1):    #X direction
        for k in range(number_of_latitudes-1): #Y direction
            v1 = vertices[l,k+1]
            v2 = vertices[l,k]
            v3 = vertices[l+1,k+1]
            v4 = vertices[l+1,k]

            #Note, this is different to the ferrit2sww code
            #since the order of the lats is reversed.
            volumes.append([v1,v3,v2]) #Upper element
            volumes.append([v4,v2,v3]) #Lower element

    volumes = array(volumes)

    geo_ref = Geo_reference(refzone,min(x),min(y))
    geo_ref.write_NetCDF(outfile)

    # This will put the geo ref in the middle
    #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)


    if verbose:
        print '------------------------------------------------'
        print 'More Statistics:'
        print '  Extent (/lon):'
        print '    x in [%f, %f], len(lat) == %d'\
              %(min(x), max(x),
                len(x))
        print '    y in [%f, %f], len(lon) == %d'\
              %(min(y), max(y),
                len(y))
        print 'geo_ref: ',geo_ref

    z = resize(bath_grid,outfile.variables['z'][:].shape)
    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
    outfile.variables['z'][:] = z
    outfile.variables['elevation'][:] = z  #FIXME HACK
    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64

    # do this to create an ok sww file.
    #outfile.variables['time'][:] = [0]   #Store time relative
    #outfile.variables['stage'] = z
    # put the next line up in the code after outfile.order = 1
    #number_of_times = 1

    stage = outfile.variables['stage']
    xmomentum = outfile.variables['xmomentum']
    ymomentum = outfile.variables['ymomentum']

    outfile.variables['time'][:] = times   #Store time relative

    if verbose: print 'Converting quantities'
    n = number_of_times
    for j in range(number_of_times):
        # load in files
        elevation_meta, elevation_grid = \
           _read_asc(elevation_dir + os.sep + elevation_files[j])

        _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
        _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])

        #print "elevation_grid",elevation_grid
        #cut matrix to desired size
        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
        #print "elevation_grid",elevation_grid
        # handle missing values
        missing = (elevation_grid == elevation_meta['NODATA_value'])
        if sometrue (missing):
            if fail_on_NaN:
                msg = 'File %s contains missing values'\
                      %(elevation_files[j])
                raise msg
            else:
                elevation_grid = elevation_grid*(missing==0) + \
                                 missing*elevation_NaN_filler


        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
        i = 0
        for k in range(number_of_latitudes):      #Y direction
            for l in range(number_of_longitudes): #X direction
                w = zscale*elevation_grid[k,l] + mean_stage
                stage[j,i] = w
                h = w - z[i]
                xmomentum[j,i] = u_momentum_grid[k,l]*h
                ymomentum[j,i] = v_momentum_grid[k,l]*h
                i += 1
    outfile.close()

def _get_min_max_indexes(latitudes,longitudes,
                        minlat=None, maxlat=None,
                        minlon=None, maxlon=None):
    """
    return max, min indexes (for slicing) of the lat and long arrays to cover the area
    specified with min/max lat/long

    Think of the latitudes and longitudes describing a 2d surface.
    The area returned is, if possible, just big enough to cover the
    inputed max/min area. (This will not be possible if the max/min area
    has a section outside of the latitudes/longitudes area.)

    assume latitudess & longitudes are sorted,
    long - from low to high (west to east, eg 148 - 151)
    lat - from high to low (north to south, eg -35 - -38)
    """

    # reverse order of lat, so it's in ascending order
    latitudes.reverse()
    largest_lat_index = len(latitudes)-1
    #Cut out a smaller extent.
    if minlat == None:
        lat_min_index = 0
    else:
        lat_min_index = searchsorted(latitudes, minlat)-1
        if lat_min_index <0:
            lat_min_index = 0


    if maxlat == None:
        lat_max_index = largest_lat_index #len(latitudes)
    else:
        lat_max_index = searchsorted(latitudes, maxlat)
        if lat_max_index > largest_lat_index:
            lat_max_index = largest_lat_index

    if minlon == None:
        lon_min_index = 0
    else:
        lon_min_index = searchsorted(longitudes, minlon)-1
        if lon_min_index <0:
            lon_min_index = 0

    if maxlon == None:
        lon_max_index = len(longitudes)
    else:
        lon_max_index = searchsorted(longitudes, maxlon)

    #Take into account that the latitude list was reversed
    latitudes.reverse() # Python passes by reference, need to swap back
    lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
                                   largest_lat_index - lat_min_index
    lat_max_index = lat_max_index + 1 # taking into account how slicing works
    lon_max_index = lon_max_index + 1 # taking into account how slicing works

    return lat_min_index, lat_max_index, lon_min_index, lon_max_index


def _read_asc(filename, verbose=False):
    """Read esri file from the following ASCII format (.asc)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    """

    datafile = open(filename)

    if verbose: print 'Reading DEM from %s' %(filename)
    lines = datafile.readlines()
    datafile.close()

    if verbose: print 'Got', len(lines), ' lines'

    ncols = int(lines.pop(0).split()[1].strip())
    nrows = int(lines.pop(0).split()[1].strip())
    xllcorner = float(lines.pop(0).split()[1].strip())
    yllcorner = float(lines.pop(0).split()[1].strip())
    cellsize = float(lines.pop(0).split()[1].strip())
    NODATA_value = float(lines.pop(0).split()[1].strip())

    assert len(lines) == nrows

    #Store data
    grid = []

    n = len(lines)
    for i, line in enumerate(lines):
        cells = line.split()
        assert len(cells) == ncols
        grid.append(array([float(x) for x in cells]))
    grid = array(grid)

    return {'xllcorner':xllcorner,
            'yllcorner':yllcorner,
            'cellsize':cellsize,
            'NODATA_value':NODATA_value}, grid