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

"""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'):
        from Scientific.IO.NetCDF import NetCDFFile
        from Numeric import Int, Float, Float32 

        self.precision = Float32 #Use single 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 = '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

                
        #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']
        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):
    """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]
    xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax]
    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]    

    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          

    #print times
    #print latitudes
    #print longitudes

    #print 'MIN', min(min(min(amplitudes)))
    #print 'MAX', max(max(max(amplitudes)))    

    #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
    
    #print file_h.dimensions.keys()
    #print file_h.variables.keys()    

    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
    #volumes = zeros(number_of_volumes, Int)
    i = 0

    #Check zone boundaries
    refzone, _, _ = redfearn(latitudes[0],longitudes[0])    

    vertices = {}
    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):
        for k in range(number_of_latitudes-1):
            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

    outfile.variables['x'][:] = x - xllcorner
    outfile.variables['y'][:] = y - yllcorner
    outfile.variables['z'][:] = 0.0
    outfile.variables['elevation'][:] = 0.0
    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):
            for l in range(number_of_longitudes):                    
                h = zscale*amplitudes[j,k,l]/100 + mean_stage
                stage[j,i] = h
                xmomentum[j,i] = xspeed[j,k,l]/100*h
                ymomentum[j,i] = yspeed[j,k,l]/100*h                
                i += 1
        
    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)]


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