source: inundation/ga/storm_surge/pyvolution/data_manager.py @ 1097

"""Functions to store and retrieve data for the Shallow Water Wave equation.
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.


Formats used within ANUGA:

.sww: Netcdf format for storing model output

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

.asc: ASCII foramt of regular DEMs as output from ArcView
.dem: NetCDF representation of regular DEM data

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

.nc: Native ferret NetCDF format

A typical dataflow can be described as follows

Manually created files:
ASC, PRJ:     Digital elevation models (gridded)
TSH:          Triangular meshes (e.g. dreated 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


"""

from Numeric import concatenate


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
        self.number_of_volumes = len(domain)
        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)
    """


    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

            #Reference point
            #Start time in seconds since the epoch (midnight 1/1/1970)
            fid.starttime = domain.starttime
            fid.xllcorner = domain.xllcorner
            fid.yllcorner = domain.yllcorner
            fid.zone = str(domain.zone) #FIXME: ?

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

            if domain.smooth is True:
                fid.createDimension('number_of_points', len(domain.vertexlist))
            else:
                fid.createDimension('number_of_points', 3*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',))

            #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[:] = 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


        #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



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

                #FIXME: Make this general (see below)
                if name == 'stage':
                    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 dem2pts(basename_in, 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 ..........

    Convert to NetCDF pts format which is

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

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

    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'

    #Georeferencing
    outfile.zone = zone
    outfile.xllcorner = xllcorner #Easting of lower left corner
    outfile.yllcorner = yllcorner #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
    outfile.createDimension('number_of_points', nrows*ncols)
    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']

    #Store data
    #FIXME: Could perhaps be faster using array operations
    for i in range(nrows):
        if verbose: print 'Processing row %d of %d' %(i, nrows)

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

            x = j*cellsize
            points[index, :] = [x,y]
            elevation[index] = dem_elevation[i, j]

    infile.close()
    outfile.close()


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
    for i, line in enumerate(lines[6:]):
        fields = line.split()
        if verbose: 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, 
               elevation = -100): #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
    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)   

    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']

    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'
    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

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


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


    if verbose: print 'Store to SWW file %s' %swwname
    # NetCDF file definition
    outfile = NetCDFFile(swwname, 'w')
        
    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'Converted from Ferret files: %s, %s, %s'\
                          %(basename_in + '_ha.nc',
                            basename_in + '_ua.nc',
                            basename_in + '_va.nc')


    #For sww compatibility
    outfile.smoothing = 'Yes' 
    outfile.order = 1
            
    #Start time in seconds since the epoch (midnight 1/1/1970)
    outfile.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', 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


    #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:
        pass
        #FIXME: z should be obtained from MOST and passed in here
        
    outfile.variables['x'][:] = x - xllcorner
    outfile.variables['y'][:] = y - yllcorner
    outfile.variables['z'][:] = z
    outfile.variables['elevation'][:] = z  #FIXME HACK
    outfile.variables['time'][:] = times
    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
    


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


    
    for j in range(len(times)):
        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
                xmomentum[j,i] = uspeed[j,k,l]/100*h
                ymomentum[j,i] = vspeed[j,k,l]/100*h
                i += 1


    if verbose:
        x = outfile.variables['x'][:]
        y = outfile.variables['y'][:]        
        
        print '------------------------------------------------'
        print 'Statistics of output file:'
        print '  Name: %s' %swwname
        print '  Reference:'
        print '    Lower left corner: [%f, %f]'\
              %(xllcorner, yllcorner)
        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']:
            q = outfile.variables[name][:].flat
            print '    %s in [%f, %f]' %(name, min(q), max(q))
            
        

        
    outfile.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,t=None,fail_if_NaN=True,NaN_filler=0):
    """Read sww Net CDF file containing Shallow Water Wave simulation

    Quantities stage, elevation, xmomentum and ymomentum.

    The momentum is not always stored.

    """
    NaN=9.969209968386869e+036
    from Scientific.IO.NetCDF import NetCDFFile
    from domain import Domain
    from Numeric import asarray, transpose
    #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    
#################################
    xllcorner = fid.xllcorner[0]
    yllcorner = fid.yllcorner[0]
    starttime = fid.starttime[0]
    zone = fid.zone
    volumes = fid.variables['volumes'][:] #Connectivity
    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
#
    conserved_quantities = []
    interpolated_quantities = {}
    other_quantities = []
#
    #print '    interpolating 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')
#
    #print other_quantities
    #print conserved_quantities
    #print '    building domain'
    domain = Domain(coordinates, volumes,\
                    conserved_quantities = conserved_quantities,\
                    other_quantities = other_quantities,zone=zone,\
                    xllcorner=xllcorner, yllcorner=yllcorner)
    domain.starttime=starttime
    domain.time=t
    for quantity in other_quantities:
        X = fid.variables[quantity][:]
        #print quantity
        #print 'max(X)'
        #print max(X)
        #print 'max(X)==NaN'
        #print max(X)==NaN
        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
        domain.set_quantity(quantity,X)
#
    for quantity in conserved_quantities:
        X = interpolated_quantities[quantity]
        #print quantity
        #print 'max(X)'
        #print max(X)
        #print 'max(X)==NaN'
        #print max(X)==NaN
        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
        domain.set_quantity(quantity,X)
    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)



#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.
    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

    # 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