Changeset 6080

Timestamp:

Dec 16, 2008, 10:15:44 AM (15 years ago)

Author:

rwilson

Message:

Fixed a few bugs, pep8 formatted.

File:

• anuga_core/source/anuga/shallow_water/data_manager.py (modified) (326 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -49 +49 @@
 """
 
-# This file was reverted from changeset:5484 to changeset:5470 on 10th July 
+# This file was reverted from changeset:5484 to changeset:5470 on 10th July
 # by Ole.
 
+import os, sys
+import csv
 import exceptions
-class TitleValueError(exceptions.Exception): pass
-class DataMissingValuesError(exceptions.Exception): pass
-class DataFileNotOpenError(exceptions.Exception): pass
-class DataTimeError(exceptions.Exception): pass
-class DataDomainError(exceptions.Exception): pass
-class NewQuantity(exceptions.Exception): pass
-
-
-
-import csv
-import os, sys
+import string
 import shutil
 from struct import unpack
 import array as p_array
-#import time, os
 from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
-
 
 from Numeric import concatenate, array, Float, Int, Int32, resize, \
@@ -76 +66 @@
      argmax, alltrue, shape, Float32, size
 
-import string
-
 from Scientific.IO.NetCDF import NetCDFFile
-#from shutil import copy
 from os.path import exists, basename, join
 from os import getcwd
-
 
 from anuga.coordinate_transforms.redfearn import redfearn, \
@@ -90 +76 @@
 from anuga.geospatial_data.geospatial_data import Geospatial_data,\
      ensure_absolute
-from anuga.config import minimum_storable_height as default_minimum_storable_height
+from anuga.config import minimum_storable_height as \
+     default_minimum_storable_height
 from anuga.config import max_float
 from anuga.utilities.numerical_tools import ensure_numeric,  mean
@@ -100 +87 @@
 from anuga.abstract_2d_finite_volumes.util import get_revision_number, \
      remove_lone_verts, sww2timeseries, get_centroid_values
-     
+
 from anuga.abstract_2d_finite_volumes.neighbour_mesh import segment_midpoints
 from anuga.load_mesh.loadASCII import export_mesh_file
@@ -106 +93 @@
 
 
+######
+# Exception classes
+######
+
+class TitleValueError(exceptions.Exception): pass
+class DataMissingValuesError(exceptions.Exception): pass
+class DataFileNotOpenError(exceptions.Exception): pass
+class DataTimeError(exceptions.Exception): pass
+class DataDomainError(exceptions.Exception): pass
+class NewQuantity(exceptions.Exception): pass
+
+
+######
 # formula mappings
+######
 
 quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
@@ -114 +115 @@
 
 
-    
+##
+# @brief Convert a possible filename into a standard form.
+# @param s Filename to process.
+# @return The new filename string.
 def make_filename(s):
     """Transform argument string into a Sexsuitable filename
@@ -128 +132 @@
 
 
+##
+# @brief Check that a specified filesystem directory path exists.
+# @param path The dirstory path to check.
+# @param verbose True if this function is to be verbose.
+# @note If directory path doesn't exist, it will be created.
 def check_dir(path, verbose=None):
     """Check that specified path exists.
@@ -141 +150 @@
     RETURN VALUE:
         Verified path including trailing separator
-
     """
 
@@ -151 +159 @@
         unix = 1
 
-
+    # add terminal separator, if it's not already there
     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 == ''):
+        path = path + os.sep
+
+    # expand ~ or ~username in path
+    path = os.path.expanduser(path)
+
+    # create directory if required
+    if not (os.access(path, os.R_OK and os.W_OK) or path == ''):
         try:
-            exitcode=os.mkdir(path)
+            exitcode = os.mkdir(path)
 
             # Change access rights if possible
-            #
             if unix:
-                exitcode=os.system('chmod 775 '+path)
+                exitcode = os.system('chmod 775 ' + path)
             else:
-                pass  # FIXME: What about acces rights under Windows?
+                pass  # FIXME: What about access rights under Windows?
 
             if verbose: print 'MESSAGE: Directory', path, 'created.'
-
         except:
             print 'WARNING: Directory', path, 'could not be created.'
@@ -174 +183 @@
                 path = '/tmp/'
             else:
-                path = 'C:'
-
-            print 'Using directory %s instead' %path
-
-    return(path)
-
-
-
+                path = 'C:' + os.sep
+
+            print "Using directory '%s' instead" % path
+
+    return path
+
+
+##
+# @brief Delete directory and all sub-directories.
+# @param path Path to the directory to delete.
 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):
@@ -199 +207 @@
                     os.remove(X)
                 except:
-                    print "Could not remove file %s" %X
+                    print "Could not remove file %s" % X
 
         os.rmdir(path)
-        
-        
-# ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007    
-def rmgeneric(path, __func__,verbose=False):
+
+
+##
+# @brief ??
+# @param path 
+# @param __func__ 
+# @param verbose True if this function is to be verbose.
+# @note ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007
+def rmgeneric(path, func, verbose=False):
     ERROR_STR= """Error removing %(path)s, %(error)s """
 
     try:
-        __func__(path)
+        func(path)
         if verbose: print 'Removed ', path
     except OSError, (errno, strerror):
         print ERROR_STR % {'path' : path, 'error': strerror }
-            
-def removeall(path,verbose=False):
-
+
+
+##
+# @brief Remove directory and all sub-directories.
+# @param path Filesystem path to directory to remove.
+# @param verbose True if this function is to be verbose.
+def removeall(path, verbose=False):
     if not os.path.isdir(path):
         return
-    
-    files=os.listdir(path)
-
-    for x in files:
-        fullpath=os.path.join(path, x)
+
+    for x in os.listdir(path):
+        fullpath = os.path.join(path, x)
         if os.path.isfile(fullpath):
-            f=os.remove
+            f = os.remove
             rmgeneric(fullpath, f)
         elif os.path.isdir(fullpath):
             removeall(fullpath)
-            f=os.rmdir
-            rmgeneric(fullpath, f,verbose)
-
-
-
+            f = os.rmdir
+            rmgeneric(fullpath, f, verbose)
+
+
+##
+# @brief Create a standard filename.
+# @param datadir Directory where file is to be created.
+# @param filename Filename 'stem'.
+# @param format Format of the file, becomes filename extension.
+# @param size Size of file, becomes part of filename.
+# @param time Time (float), becomes part of filename.
+# @return The complete filename path, including directory.
+# @note The containing directory is created, if necessary.
 def create_filename(datadir, filename, format, size=None, time=None):
-
-    import os
-    #from anuga.config import data_dir
-
     FN = check_dir(datadir) + filename
 
     if size is not None:
-        FN += '_size%d' %size
+        FN += '_size%d' % size
 
     if time is not None:
-        FN += '_time%.2f' %time
+        FN += '_time%.2f' % time
 
     FN += '.' + format
+
     return FN
 
 
+##
+# @brief Get all files with a standard name and a given set of attributes.
+# @param datadir Directory files must be in.
+# @param filename Filename stem.
+# @param format Filename extension.
+# @param size Filename size.
+# @return A list of fielnames (including directory) that match the attributes.
 def get_files(datadir, filename, format, size):
     """Get all file (names) with given name, size and format
@@ -256 +283 @@
     import glob
 
-    import os
-    #from anuga.config import data_dir
-
     dir = check_dir(datadir)
-
-    pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format)
+    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
+##
+# @brief 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)'''
+    ##
+    # @brief Instantiate this instance.
+    # @param domain 
+    # @param extension 
+    # @param mode The mode of the underlying file.
+    def __init__(self, domain, extension, mode='w'):
+        assert mode in ['r', 'w', 'a'], \
+               "Mode %s must be either:\n" % mode + \
+               "   'w' (write)\n" + \
+               "   'r' (read)\n" + \
+               "   'a' (append)"
 
         #Create filename
@@ -282 +311 @@
                                         domain.get_name(), extension)
 
-        #print 'F', self.filename
         self.timestep = 0
         self.domain = domain
-        
-
 
         # Exclude ghosts in case this is a parallel domain
-        self.number_of_nodes = domain.number_of_full_nodes        
+        self.number_of_nodes = domain.number_of_full_nodes
         self.number_of_volumes = domain.number_of_full_triangles
-        #self.number_of_volumes = len(domain)        
-
-
-
+        #self.number_of_volumes = len(domain)
 
         #FIXME: Should we have a general set_precision function?
 
 
-
-#Class for storing output to e.g. visualisation
+##
+# @brief Class for storing output to e.g. visualisation
 class Data_format_sww(Data_format):
     """Interface to native NetCDF format (.sww) for storing model output
@@ -312 +335 @@
     """
 
-
-    def __init__(self, domain, mode = 'w',\
-                 max_size = 2000000000,
-                 recursion = False):
+    ##
+    # @brief Instantiate this instance.
+    # @param domain ??
+    # @param mode Mode of the underlying data file.
+    # @param max_size ??
+    # @param recursion ??
+    # @note Prepare the undelying data file if mode is 'w'.
+    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 for quantities
+        self.recursion = recursion
+        self.mode = mode
         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)
-
         if hasattr(domain, 'minimum_storable_height'):
             self.minimum_storable_height = domain.minimum_storable_height
@@ -334 +358 @@
             self.minimum_storable_height = default_minimum_storable_height
 
+        # call owning constructor
+        Data_format.__init__(self, domain, 'sww', mode)
+
         # NetCDF file definition
         fid = NetCDFFile(self.filename, mode)
-
         if mode == 'w':
-            description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting'
+            description = 'Output from anuga.abstract_2d_finite_volumes ' \
+                          'suitable for plotting'
             self.writer = Write_sww()
             self.writer.store_header(fid,
@@ -355 +382 @@
             if domain.quantities_to_be_monitored is not None:
                 fid.createDimension('singleton', 1)
-                fid.createDimension('two', 2)                
+                fid.createDimension('two', 2)
 
                 poly = domain.monitor_polygon
@@ -363 +390 @@
                     fid.createVariable('extrema.polygon',
                                        self.precision,
-                                       ('polygon_length',
-                                        'two'))
-                    fid.variables['extrema.polygon'][:] = poly                                    
-
-                    
+                                       ('polygon_length', 'two'))
+                    fid.variables['extrema.polygon'][:] = poly
+
                 interval = domain.monitor_time_interval
                 if interval is not None:
@@ -375 +400 @@
                     fid.variables['extrema.time_interval'][:] = interval
 
-                
                 for q in domain.quantities_to_be_monitored:
-                    #print 'doing', q
-                    fid.createVariable(q+'.extrema', self.precision,
+                    fid.createVariable(q + '.extrema', self.precision,
                                        ('numbers_in_range',))
-                    fid.createVariable(q+'.min_location', self.precision,
+                    fid.createVariable(q + '.min_location', self.precision,
                                        ('numbers_in_range',))
-                    fid.createVariable(q+'.max_location', self.precision,
+                    fid.createVariable(q + '.max_location', self.precision,
                                        ('numbers_in_range',))
-                    fid.createVariable(q+'.min_time', self.precision,
+                    fid.createVariable(q + '.min_time', self.precision,
                                        ('singleton',))
-                    fid.createVariable(q+'.max_time', self.precision,
+                    fid.createVariable(q + '.max_time', self.precision,
                                        ('singleton',))
 
-                    
         fid.close()
 
-
+    ##
+    # @brief Store connectivity data into the underlying data file.
     def store_connectivity(self):
         """Specialisation of store_connectivity for net CDF format
@@ -401 +424 @@
 
         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
+        # append to the NetCDF file
+        fid = NetCDFFile(self.filename, 'a')
 
         # Get the variables
@@ -413 +435 @@
         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,Y,Z,V = Q.get_vertex_values(xy=True, precision=self.precision)
+
+        # store the connectivity data
         points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
         self.writer.store_triangulation(fid,
@@ -427 +447 @@
                                         V.astype(volumes.typecode()),
                                         Z,
-                                        points_georeference= \
-                                        domain.geo_reference)
-
-        # Close
+                                        points_georeference=\
+                                            domain.geo_reference)
+
         fid.close()
 
-
+    ##
+    # @brief Store time and named quantities to the underlying data file.
+    # @param names The names of the quantities to store.
+    # @note If 'names' not supplied, store a standard set of quantities.
     def store_timestep(self, names=None):
         """Store time and named quantities to file
         """
-        
+
         from Scientific.IO.NetCDF import NetCDFFile
         import types
         from time import sleep
         from os import stat
-
         from Numeric import choose
-
 
         if names is None:
             # Standard shallow water wave equation quantitites in ANUGA
             names = ['stage', 'xmomentum', 'ymomentum']
-        
-        # Get NetCDF       
+
+        # Get NetCDF
         retries = 0
         file_open = False
@@ -459 +479 @@
                 # This could happen if someone was reading the file.
                 # In that case, wait a while and try again
-                msg = 'Warning (store_timestep): File %s could not be opened'\
-                      %self.filename
-                msg += ' - trying step %s again' %self.domain.time
+                msg = 'Warning (store_timestep): File %s could not be opened' \
+                      % self.filename
+                msg += ' - trying step %s again' % self.domain.time
                 print msg
                 retries += 1
@@ -469 +489 @@
 
         if not file_open:
-            msg = 'File %s could not be opened for append' %self.filename
+            msg = 'File %s could not be opened for append' % self.filename
             raise DataFileNotOpenError, msg
-
-
 
         # Check to see if the file is already too big:
         time = fid.variables['time']
-        i = len(time)+1
+        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):
+        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 domain.filename and domain.starttime.
@@ -487 +505 @@
             # Write a filename addon that won't break swollens reader
             # (10.sww is bad)
-            filename_ext = '_time_%s'%self.domain.time
+            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.get_name()
             if not self.recursion:
-                self.domain.set_name(old_domain_filename+filename_ext)
-
+                self.domain.set_name(old_domain_filename + filename_ext)
 
             # Change the domain starttime to the current time
@@ -502 +519 @@
 
             # 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)
+            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
+                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
@@ -520 +537 @@
             #restore the old starttime and filename
             self.domain.starttime = old_domain_starttime
-            self.domain.set_name(old_domain_filename)            
+            self.domain.set_name(old_domain_filename)
         else:
             self.recursion = False
@@ -534 +551 @@
                 names = [names]
 
-            if 'stage' in names and 'xmomentum' in names and \
-               'ymomentum' in names:
-
+            if 'stage' in names \
+               and 'xmomentum' in names \
+               and 'ymomentum' in names:
                 # Get stage, elevation, depth and select only those
                 # values where minimum_storable_height is exceeded
                 Q = domain.quantities['stage']
-                A, _ = Q.get_vertex_values(xy = False,
-                                           precision = self.precision)
+                A, _ = Q.get_vertex_values(xy=False, precision=self.precision)
                 z = fid.variables['elevation']
 
-                storable_indices = A-z[:] >= self.minimum_storable_height
+                storable_indices = (A-z[:] >= self.minimum_storable_height)
                 stage = choose(storable_indices, (z[:], A))
-                
+
                 # Define a zero vector of same size and type as A
                 # for use with momenta
                 null = zeros(size(A), A.typecode())
-                
+
                 # Get xmomentum where depth exceeds minimum_storable_height
                 Q = domain.quantities['xmomentum']
-                xmom, _ = Q.get_vertex_values(xy = False,
-                                              precision = self.precision)
+                xmom, _ = Q.get_vertex_values(xy=False,
+                                              precision=self.precision)
                 xmomentum = choose(storable_indices, (null, xmom))
-                
+
 
                 # Get ymomentum where depth exceeds minimum_storable_height
                 Q = domain.quantities['ymomentum']
-                ymom, _ = Q.get_vertex_values(xy = False,
-                                              precision = self.precision)
-                ymomentum = choose(storable_indices, (null, ymom))                
-                
-                # Write quantities to NetCDF
-                self.writer.store_quantities(fid, 
+                ymom, _ = Q.get_vertex_values(xy=False,
+                                              precision=self.precision)
+                ymomentum = choose(storable_indices, (null, ymom))
+
+                # Write quantities to underlying data  file
+                self.writer.store_quantities(fid,
                                              time=self.domain.time,
                                              sww_precision=self.precision,
@@ -572 +588 @@
                                              ymomentum=ymomentum)
             else:
-                msg = 'Quantities stored must be: stage, xmomentum, ymomentum.'
-                msg += ' Instead I got: ' + str(names)
+                msg = 'Quantities stored must be: stage, xmomentum, ymomentum, '
+                msg += 'but I got: ' + str(names)
                 raise Exception, msg
-            
-
 
             # Update extrema if requested
@@ -582 +596 @@
             if domain.quantities_to_be_monitored is not None:
                 for q, info in domain.quantities_to_be_monitored.items():
-
                     if info['min'] is not None:
                         fid.variables[q + '.extrema'][0] = info['min']
-                        fid.variables[q + '.min_location'][:] =\
+                        fid.variables[q + '.min_location'][:] = \
                                         info['min_location']
                         fid.variables[q + '.min_time'][0] = info['min_time']
-                        
+
                     if info['max'] is not None:
                         fid.variables[q + '.extrema'][1] = info['max']
-                        fid.variables[q + '.max_location'][:] =\
+                        fid.variables[q + '.max_location'][:] = \
                                         info['max_location']
                         fid.variables[q + '.max_time'][0] = info['max_time']
-
-            
 
             # Flush and close
@@ -602 +613 @@
 
 
-
-# Class for handling checkpoints data
+##
+# @brief Class for handling checkpoints data
 class Data_format_cpt(Data_format):
     """Interface to native NetCDF format (.cpt)
     """
 
-
-    def __init__(self, domain, mode = 'w'):
+    ##
+    # @brief Initialize this instantiation.
+    # @param domain ??
+    # @param mode Mode of underlying data file (default WRITE).
+    def __init__(self, domain, mode='w'):
         from Scientific.IO.NetCDF import NetCDFFile
         from Numeric import Int, Float, Float
@@ -617 +631 @@
         Data_format.__init__(self, domain, 'sww', mode)
 
-
         # NetCDF file definition
         fid = NetCDFFile(self.filename, mode)
-
         if mode == 'w':
             #Create new file
@@ -646 +658 @@
                                                 'number_of_vertices'))
 
-            fid.createVariable('time', self.precision,
-                               ('number_of_timesteps',))
+            fid.createVariable('time', self.precision, ('number_of_timesteps',))
 
             #Allocate space for all quantities
@@ -655 +666 @@
                                     'number_of_points'))
 
-        #Close
         fid.close()
 
-
+    ##
+    # @brief Store connectivity data to underlying data file.
     def store_checkpoint(self):
-        """
-        Write x,y coordinates of triangles.
+        """Write x,y coordinates of triangles.
         Write connectivity (
         constituting
@@ -668 +678 @@
 
         from Scientific.IO.NetCDF import NetCDFFile
-
         from Numeric import concatenate
 
@@ -674 +683 @@
 
         #Get NetCDF
-        fid = NetCDFFile(self.filename, 'a')  #Open existing file for append
+        fid = NetCDFFile(self.filename, 'a')
 
         # Get the variables
@@ -684 +693 @@
         # 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,Y,Z,V = Q.get_vertex_values(xy=True, precision=self.precision)
 
         x[:] = X.astype(self.precision)
@@ -695 +701 @@
         volumes[:] = V
 
-        #Close
         fid.close()
 
-
+    ##
+    # @brief Store tiem and named quantities to underlying data file.
+    # @param name 
     def store_timestep(self, name):
         """Store time and named quantity to file
         """
+
         from Scientific.IO.NetCDF import NetCDFFile
         from time import sleep
@@ -710 +718 @@
         while not file_open and retries < 10:
             try:
-                fid = NetCDFFile(self.filename, 'a') #Open existing file
+                fid = NetCDFFile(self.filename, 'a')
             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'
+                msg = 'Warning (store_timestep): File %s could not be opened' \
+                      ' - trying again' % self.filename
                 print msg
                 retries += 1
@@ -724 +731 @@
 
         if not file_open:
-            msg = 'File %s could not be opened for append' %self.filename
-            raise DataFileNotOPenError, msg
-
+            msg = 'File %s could not be opened for append' % self.filename
+            raise DataFileNotOpenError, msg
 
         domain = self.domain
@@ -740 +746 @@
         # Get quantity
         Q = domain.quantities[name]
-        A,V = Q.get_vertex_values(xy=False,
-                                  precision = self.precision)
+        A,V = Q.get_vertex_values(xy=False, precision=self.precision)
 
         stage[i,:] = A.astype(self.precision)
@@ -750 +755 @@
 
 
-#### NED is national exposure database (name changed to NEXIS)
-    
+##
+# @brief Class for National Exposure Database storage (NEXIS).
+
 LAT_TITLE = 'LATITUDE'
 LONG_TITLE = 'LONGITUDE'
 X_TITLE = 'x'
 Y_TITLE = 'y'
+
 class Exposure_csv:
+
+    ##
+    # @brief Instantiate this instance.
+    # @param file_name Name of underlying data file.
+    # @param latitude_title ??
+    # @param longitude_title ??
+    # @param is_x_y_locations ??
+    # @param x_title ??
+    # @param y_title ??
+    # @param refine_polygon ??
+    # @param title_check_list ??
     def __init__(self,file_name, latitude_title=LAT_TITLE,
                  longitude_title=LONG_TITLE, is_x_y_locations=None,
@@ -772 +790 @@
            each column is a 'set' of data
 
-        Feel free to use/expand it to read other csv files. 
-           
-           
+        Feel free to use/expand it to read other csv files.
+
         It is not for adding and deleting rows
-        
+
         Can geospatial handle string attributes? It's not made for them.
         Currently it can't load and save string att's.
@@ -785 +802 @@
 
         The location info is in the geospatial attribute.
-        
-        
         """
+
         self._file_name = file_name
         self._geospatial = None #
@@ -794 +810 @@
         #The keys are the column titles.
         #The values are lists of column data
-        
+
         # self._title_index_dic is a dictionary.
         #The keys are the column titles.
@@ -801 +817 @@
             csv2dict(self._file_name, title_check_list=title_check_list)
         try:
-            #Have code here that handles caps or lower 
+            #Have code here that handles caps or lower
             lats = self._attribute_dic[latitude_title]
             longs = self._attribute_dic[longitude_title]
-            
         except KeyError:
             # maybe a warning..
@@ -812 +827 @@
             pass
         else:
-            self._geospatial = Geospatial_data(latitudes = lats,
-                 longitudes = longs)
+            self._geospatial = Geospatial_data(latitudes=lats,
+                                               longitudes=longs)
 
         if is_x_y_locations is True:
@@ -828 +843 @@
             else:
                 self._geospatial = Geospatial_data(data_points=points)
-            
+
         # create a list of points that are in the refining_polygon
         # described by a list of indexes representing the points
 
+    ##
+    # @brief Create a comparison method.
+    # @param self This object.
+    # @param other The other object.
+    # @return True if objects are 'same'.
     def __cmp__(self, other):
-        #print "self._attribute_dic",self._attribute_dic
-        #print "other._attribute_dic",other._attribute_dic
-        #print "self._title_index_dic", self._title_index_dic
-        #print "other._title_index_dic", other._title_index_dic
-        
-        #check that a is an instance of this class
+        #check that 'other' is an instance of this class
         if isinstance(self, type(other)):
             result = cmp(self._attribute_dic, other._attribute_dic)
-            if result <>0:
+            if result <> 0:
                 return result
-            # The order of the columns is important. Therefore.. 
+
+            # The order of the columns is important. Therefore..
             result = cmp(self._title_index_dic, other._title_index_dic)
-            if result <>0:
+            if result <> 0:
                 return result
-            for self_ls, other_ls in map(None,self._attribute_dic, \
-                    other._attribute_dic):
+            for self_ls, other_ls in map(None, self._attribute_dic,
+                                         other._attribute_dic):
                 result = cmp(self._attribute_dic[self_ls],
                              other._attribute_dic[other_ls])
-                if result <>0:
+                if result <> 0:
                     return result
             return 0
         else:
             return 1
-    
-
+
+    ##
+    # @brief Get a list of column values given a column name.
+    # @param column_name The name of the column to get values from.
+    # @param use_refind_polygon Unused??
     def get_column(self, column_name, use_refind_polygon=False):
         """
@@ -865 +884 @@
         do this to change a list of strings to a list of floats
         time = [float(x) for x in time]
-        
+
         Not implemented:
         if use_refind_polygon is True, only return values in the
         refined polygon
         """
+
         if not self._attribute_dic.has_key(column_name):
-            msg = 'Therer is  no column called %s!' %column_name
+            msg = 'There is no column called %s!' % column_name
             raise TitleValueError, msg
+
         return self._attribute_dic[column_name]
 
-
-    def get_value(self, value_column_name,
-                  known_column_name,
-                  known_values,
-                  use_refind_polygon=False):
+    ##
+    # @brief ??
+    # @param value_column_name ??
+    # @param known_column_name ??
+    # @param known_values ??
+    # @param use_refind_polygon ??
+    def get_value(self, value_column_name, known_column_name,
+                  known_values, use_refind_polygon=False):
         """
         Do linear interpolation on the known_colum, using the known_value,
         to return a value of the column_value_name.
         """
+
         pass
 
-
+    ##
+    # @brief Get a geospatial object that describes the locations.
+    # @param use_refind_polygon Unused??
     def get_location(self, use_refind_polygon=False):
         """
@@ -893 +920 @@
 
         Note, if there is not location info, this returns None.
-        
+
         Not implemented:
         if use_refind_polygon is True, only return values in the
         refined polygon
         """
+
         return self._geospatial
 
+    ##
+    # @brief Add column to 'end' of CSV data.
+    # @param column_name The new column name.
+    # @param column_values The new column values.
+    # @param overwrite If True, overwrites last column, doesn't add at end.
     def set_column(self, column_name, column_values, overwrite=False):
         """
@@ -906 +939 @@
 
         Set overwrite to True if you want to overwrite a column.
-        
+
         Note, in column_name white space is removed and case is not checked.
         Precondition
         The column_name and column_values cannot have comma's in it.
         """
-        
+
+        # sanity checks
         value_row_count = \
-            len(self._attribute_dic[self._title_index_dic.keys()[0]])
-        if len(column_values) <> value_row_count: 
+                len(self._attribute_dic[self._title_index_dic.keys()[0]])
+        if len(column_values) <> value_row_count:
             msg = 'The number of column values must equal the number of rows.'
             raise DataMissingValuesError, msg
-        
+
+        # check new column name isn't already used, and we aren't overwriting
         if self._attribute_dic.has_key(column_name):
             if not overwrite:
-                msg = 'Column name %s already in use!' %column_name
+                msg = 'Column name %s already in use!' % column_name
                 raise TitleValueError, msg
         else:
             # New title.  Add it to the title index.
             self._title_index_dic[column_name] = len(self._title_index_dic)
+
         self._attribute_dic[column_name] = column_values
-        #print "self._title_index_dic[column_name]",self._title_index_dic 
-
+
+    ##
+    # @brief Save the exposure CSV  file.
+    # @param file_name If supplied, use this filename, not original.
     def save(self, file_name=None):
         """
         Save the exposure csv file
         """
+
         if file_name is None:
             file_name = self._file_name
-        
-        fd = open(file_name,'wb')
+
+        fd = open(file_name, 'wb')
         writer = csv.writer(fd)
-        
+
         #Write the title to a cvs file
-        line = [None]* len(self._title_index_dic)
+        line = [None] * len(self._title_index_dic)
         for title in self._title_index_dic.iterkeys():
-            line[self._title_index_dic[title]]= title
+            line[self._title_index_dic[title]] = title
         writer.writerow(line)
-        
+
         # Write the values to a cvs file
         value_row_count = \
-            len(self._attribute_dic[self._title_index_dic.keys()[0]])
+                len(self._attribute_dic[self._title_index_dic.keys()[0]])
         for row_i in range(value_row_count):
-            line = [None]* len(self._title_index_dic)
+            line = [None] * len(self._title_index_dic)
             for title in self._title_index_dic.iterkeys():
-                line[self._title_index_dic[title]]= \
+                line[self._title_index_dic[title]] = \
                      self._attribute_dic[title][row_i]
             writer.writerow(line)
 
 
+##
+# @brief Convert CSV file to a dictionary of arrays.
+# @param file_name The path to the file to read.
 def csv2array(file_name):
-    """Convert CSV files of the form
-    
+    """
+    Convert CSV files of the form:
+
     time, discharge, velocity
     0.0,  1.2,       0.0
     0.1,  3.2,       1.1
     ...
-    
+
     to a dictionary of numeric arrays.
-    
-    
+
+
     See underlying function csv2dict for more details.
-    
-    """
-    
-    
+    """
+
     X, _ = csv2dict(file_name)
-    
+
     Y = {}
     for key in X.keys():
         Y[key] = array([float(x) for x in X[key]])
-        
-    return Y    
-    
-            
+
+    return Y
+
+
+##
+# @brief Read a CSV file and convert to a dictionary of {key: column}.
+# @param file_name The path to the file to read.
+# @param title_check_list List of titles that *must* be columns in the file.
+# @return Two dicts: ({key:column}, {title:index}).
+# @note WARNING: Values are returned as strings.
 def csv2dict(file_name, title_check_list=None):
     """
-    Load in the csv as a dic, title as key and column info as value, .
+    Load in the csv as a dic, title as key and column info as value.
     Also, create a dic, title as key and column index as value,
-    to keep track of the column order. 
+    to keep track of the column order.
 
     Two dictionaries are returned.
-    
+
     WARNING: Values are returned as strings.
-    do this to change a list of strings to a list of floats
-        time = [float(x) for x in time]
-
-        
-    """
-    
-    #
+             Do this to change a list of strings to a list of floats
+                 time = [float(x) for x in time]
+    """
+
     attribute_dic = {}
     title_index_dic = {}
     titles_stripped = [] # list of titles
+
     reader = csv.reader(file(file_name))
 
@@ -1006 +1051 @@
         titles_stripped.append(title.strip())
         title_index_dic[title.strip()] = i
-    title_count = len(titles_stripped)       
-    #print "title_index_dic",title_index_dic
+    title_count = len(titles_stripped)
+
+    # check required columns
     if title_check_list is not None:
         for title_check in title_check_list:
-            #msg = "Reading error.  This row is not present ", title_check 
+            #msg = "Reading error. This row is not present ", title_check
             #assert title_index_dic.has_key(title_check), msg
             if not title_index_dic.has_key(title_check):
                 #reader.close()
-                msg = "Reading error.  This row is not present ", \
-                      title_check                     
+                msg = "Reading error. This row is not present ", title_check
                 raise IOError, msg
-                
-    
-    
-    #create a dic of colum values, indexed by column title
+
+    #create a dic of column values, indexed by column title
     for line in reader:
         if len(line) <> title_count:
             raise IOError #FIXME make this nicer
         for i, value in enumerate(line):
-            attribute_dic.setdefault(titles_stripped[i],[]).append(value)
-        
+            attribute_dic.setdefault(titles_stripped[i], []).append(value)
+
     return attribute_dic, title_index_dic
 
 
-#Auxiliary
-def write_obj(filename,x,y,z):
-    """Store x,y,z vectors into filename (obj format)
+##
+# @brief 
+# @param filename 
+# @param x 
+# @param y 
+# @param z 
+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.
@@ -1040 +1089 @@
 
        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
@@ -1052 +1099 @@
         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
+    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]))
+            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))
+        base = i * N
+        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
 
     outfile.close()
@@ -1074 +1120 @@
 ########################################################
 
+##
+# @brief Convert SWW data to OBJ data.
+# @param basefilename Stem of filename, needs size and extension added.
+# @param size The number of lines to write.
 def sww2obj(basefilename, size):
     """Convert netcdf based data output to obj
     """
+
     from Scientific.IO.NetCDF import NetCDFFile
-
     from Numeric import Float, zeros
 
-    #Get NetCDF
+    # Get NetCDF
     FN = create_filename('.', basefilename, 'sww', size)
     print 'Reading from ', FN
     fid = NetCDFFile(FN, 'r')  #Open existing file for read
-
 
     # Get the variables
@@ -1105 +1154 @@
             zz[i,j] = z[i+j*M]
 
-    #Write obj for bathymetry
+    # 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
-
+    # 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]
@@ -1121 +1168 @@
                 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)
-
-
+        write_obj(FN, xx, yy, zz)
+
+
+##
+# @brief 
+# @param basefilename Stem of filename, needs size and extension added.
 def dat2obj(basefilename):
     """Convert line based data output to obj
@@ -1135 +1184 @@
     import glob, os
     from anuga.config import data_dir
-
-
-    #Get bathymetry and x,y's
+    from Numeric import zeros, Float
+
+    # 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
@@ -1147 +1194 @@
     z = zeros((M,3), Float)
 
-    ##i = 0
     for i, line in enumerate(lines):
         tokens = line.split()
-        values = map(float,tokens)
+        values = map(float, tokens)
 
         for j in range(3):
@@ -1157 +1203 @@
             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')
-
+    # 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()
+        lines = open(data_dir + os.sep + filename, 'r').readlines()
         assert len(lines) == M
         root, ext = os.path.splitext(filename)
 
-        #Get time from filename
+        # Get time from filename
         i0 = filename.find('_time=')
         if i0 == -1:
@@ -1191 +1231 @@
             raise DataTimeError, 'Hmmmm'
 
-
-
-        ##i = 0
         for i, line in enumerate(lines):
             tokens = line.split()
@@ -1201 +1238 @@
                 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):
+        # Write obj for variable data
+        write_obj(data_dir + os.sep + basefilename + '_time=%.4f' % t, x, y, z)
+
+
+##
+# @brief Filter data file, selecting timesteps first:step:last.
+# @param filename1 Data file to filter.
+# @param filename2 File to write filtered timesteps to.
+# @param first First timestep.
+# @param last Last timestep.
+# @param step Timestep stride.
+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
+    # Get NetCDF
     infile = NetCDFFile(filename1, 'r')  #Open existing file for read
     outfile = NetCDFFile(filename2, 'w')  #Open new file
 
-
-    #Copy dimensions
+    # Copy dimensions
     for d in infile.dimensions:
         outfile.createDimension(d, infile.dimensions[d])
 
+    # Copy variable definitions
     for name in infile.variables:
         var = infile.variables[name]
         outfile.createVariable(name, var.typecode(), var.dimensions)
 
-
-    #Copy the static variables
+    # 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
+    # Copy selected timesteps
     time = infile.variables['time']
     stage = infile.variables['stage']
@@ -1247 +1288 @@
     selection = range(first, last, step)
     for i, j in enumerate(selection):
-        print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j])
+        print 'Copying timestep %d of %d (%f)' % (j, last-first, time[j])
         newtime[i] = time[j]
         newstage[i,:] = stage[j,:]
 
-    #Close
+    # Close
     infile.close()
     outfile.close()
 
 
+##
+# @brief Return instance of class of given format using filename.
+# @param domain Data domain (eg, 'sww', etc).
+# @param mode The mode to open domain in.
+# @return A class instance of required domain and mode.
 #Get data objects
 def get_dataobject(domain, mode='w'):
@@ -1261 +1307 @@
     """
 
-    cls = eval('Data_format_%s' %domain.format)
+    cls = eval('Data_format_%s' % domain.format)
     return cls(domain, mode)
 
 
-
-
+##
+# @brief Convert DEM data  to PTS data.
+# @param basename_in Stem of input filename.
+# @param basename_out Stem of output filename.
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
+# @param use_cache 
+# @param verbose 
+# @return 
 def dem2pts(basename_in, basename_out=None,
             easting_min=None, easting_max=None,
@@ -1289 +1344 @@
     """
 
-
-
     kwargs = {'basename_out': basename_out,
               'easting_min': easting_min,
@@ -1310 +1363 @@
 
 
+##
+# @brief 
+# @param basename_in 
+# @param basename_out 
+# @param verbose 
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
 def _dem2pts(basename_in, basename_out=None, verbose=False,
             easting_min=None, easting_max=None,
@@ -1327 +1389 @@
 
     # Get NetCDF
-    infile = NetCDFFile(root + '.dem', 'r')  # Open existing netcdf file for read
+    infile = NetCDFFile(root + '.dem', 'r') 
 
     if verbose: print 'Reading DEM from %s' %(root + '.dem')
@@ -1348 +1410 @@
     units = infile.units
 
-
     # Get output file
     if basename_out == None:
@@ -1356 +1417 @@
 
     if verbose: print 'Store to NetCDF file %s' %ptsname
+
     # NetCDF file definition
     outfile = NetCDFFile(ptsname, 'w')
@@ -1361 +1423 @@
     # Create new file
     outfile.institution = 'Geoscience Australia'
-    outfile.description = 'NetCDF pts format for compact and portable storage ' +\
-                      'of spatial point data'
+    outfile.description = 'NetCDF pts format for compact and portable ' \
+                          'storage of spatial point data'
+
     # Assign default values
     if easting_min is None: easting_min = xllcorner
@@ -1383 +1446 @@
     outfile.datum = datum
     outfile.units = units
-
 
     # Grid info (FIXME: probably not going to be used, but heck)
@@ -1404 +1466 @@
         for j in range(ncols):
             x = j*cellsize + xllcorner
-            if easting_min <= x <= easting_max and \
-               northing_min <= y <= northing_max:
+            if easting_min <= x <= easting_max \
+               and northing_min <= y <= northing_max:
                 thisj = j
                 thisi = i
-                if dem_elevation_r[i,j] == NODATA_value: nn += 1
+                if dem_elevation_r[i,j] == NODATA_value:
+                    nn += 1
 
                 if k == 0:
                     i1_0 = i
                     j1_0 = j
+
                 k += 1
 
@@ -1445 +1509 @@
 
     lenv = index2-index1+1
+
     # Store data
     global_index = 0
     # for i in range(nrows):
-    for i in range(i1_0,thisi+1,1):
-        if verbose and i%((nrows+10)/10)==0:
-            print 'Processing row %d of %d' %(i, nrows)
+    for i in range(i1_0, thisi+1, 1):
+        if verbose and i % ((nrows+10)/10) == 0:
+            print 'Processing row %d of %d' % (i, nrows)
 
         lower_index = global_index
@@ -1457 +1522 @@
         no_NODATA = sum(v == NODATA_value)
         if no_NODATA > 0:
-            newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA
+            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
         else:
-            newcols = lenv # ncols_in_bounding_box
+            newcols = lenv              # ncols_in_bounding_box
 
         telev = zeros(newcols, Float)
@@ -1469 +1534 @@
         #for j in range(ncols):
         for j in range(j1_0,index2+1,1):
-
             x = j*cellsize + xllcorner
-            if easting_min <= x <= easting_max and \
-               northing_min <= y <= northing_max and \
-               dem_elevation_r[i,j] <> NODATA_value:
-                tpoints[local_index, :] = [x-easting_min,y-northing_min]
+            if easting_min <= x <= easting_max \
+               and northing_min <= y <= northing_max \
+               and dem_elevation_r[i,j] <> NODATA_value:
+                tpoints[local_index, :] = [x-easting_min, y-northing_min]
                 telev[local_index] = dem_elevation_r[i, j]
                 global_index += 1
@@ -1491 +1555 @@
 
 
-
+##
+# @brief  Return block of surface lines for each cross section
+# @param lines Iterble  of text lines to process.
+# @note BROKEN?  UNUSED?
 def _read_hecras_cross_sections(lines):
     """Return block of surface lines for each cross section
@@ -1502 +1569 @@
     reading_surface = False
     for i, line in enumerate(lines):
-
         if len(line.strip()) == 0:    #Ignore blanks
             continue
@@ -1523 +1589 @@
 
 
-
-
+##
+# @brief Convert HECRAS (.sdf) file to PTS file.
+# @param basename_in Sterm of input filename.
+# @param basename_out Sterm of output filename.
+# @param verbose True if this function is to be verbose.
 def hecras_cross_sections2pts(basename_in,
                               basename_out=None,
@@ -1531 +1600 @@
 
     Example:
-
 
 # RAS export file created on Mon 15Aug2005 11:42
@@ -1550 +1618 @@
 END HEADER:
 
-
 Only the SURFACE LINE data of the following form will be utilised
-
   CROSS-SECTION:
     STREAM ID:Southern-Wungong
@@ -1587 +1653 @@
     root = basename_in
 
-    #Get ASCII file
-    infile = open(root + '.sdf', 'r')  #Open SDF file for read
+    # Get ASCII file
+    infile = open(root + '.sdf', 'r')
 
     if verbose: print 'Reading DEM from %s' %(root + '.sdf')
@@ -1597 +1663 @@
     if verbose: print 'Converting to pts format'
 
+    # Scan through the header, picking up stuff we need.
     i = 0
     while lines[i].strip() == '' or lines[i].strip().startswith('#'):
@@ -1642 +1709 @@
     i += 1
 
-
-    #Now read all points
+    # Now read all points
     points = []
     elevation = []
@@ -1651 +1717 @@
             elevation.append(entry[2])
 
-
     msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
           %(j+1, number_of_cross_sections)
     assert j+1 == number_of_cross_sections, msg
 
-    #Get output file
+    # Get output file, write PTS data
     if basename_out == None:
         ptsname = root + '.pts'
@@ -1667 +1732 @@
     geo.export_points_file(ptsname)
 
-def export_grid(basename_in, extra_name_out = None,
-                quantities = None, # defaults to elevation
-                timestep = None,
-                reduction = None,
-                cellsize = 10,
-                number_of_decimal_places = None,
-                NODATA_value = -9999,
-                easting_min = None,
-                easting_max = None,
-                northing_min = None,
-                northing_max = None,
-                verbose = False,
-                origin = None,
-                datum = 'WGS84',
-                format = 'ers'):
-    """
-    
-    Wrapper for sww2dem. - see sww2dem to find out what most of the
-    parameters do.
+
+##
+# @brief 
+# @param basename_in 
+# @param extra_name_out 
+# @param quantities 
+# @param timestep 
+# @param reduction 
+# @param cellsize 
+# @param number_of_decimal_places 
+# @param NODATA_value 
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
+# @param verbose 
+# @param origin 
+# @param datum 
+# @param format 
+# @return 
+def export_grid(basename_in, extra_name_out=None,
+                quantities=None, # defaults to elevation
+                timestep=None,
+                reduction=None,
+                cellsize=10,
+                number_of_decimal_places=None,
+                NODATA_value=-9999,
+                easting_min=None,
+                easting_max=None,
+                northing_min=None,
+                northing_max=None,
+                verbose=False,
+                origin=None,
+                datum='WGS84',
+                format='ers'):
+    """Wrapper for sww2dem.
+    See sww2dem to find out what most of the parameters do.
 
     Quantities is a list of quantities.  Each quantity will be
@@ -1695 +1778 @@
     This function returns the names of the files produced.
 
-    It will also produce as many output files as there are input sww files. 
-    """
-    
+    It will also produce as many output files as there are input sww files.
+    """
+
     if quantities is None:
         quantities = ['elevation']
-        
+
     if type(quantities) is str:
             quantities = [quantities]
@@ -1706 +1789 @@
     # How many sww files are there?
     dir, base = os.path.split(basename_in)
-    #print "basename_in",basename_in
-    #print "base",base
-
-    iterate_over = get_all_swwfiles(dir,base,verbose)
-    
+
+    iterate_over = get_all_swwfiles(dir, base, verbose)
+
     if dir == "":
         dir = "." # Unix compatibility
-    
+
     files_out = []
-    #print 'sww_file',iterate_over
     for sww_file in iterate_over:
         for quantity in quantities:
@@ -1721 +1801 @@
                 basename_out = sww_file + '_' + quantity
             else:
-                basename_out = sww_file + '_' + quantity + '_' \
-                               + extra_name_out
-#            print "basename_out", basename_out
-        
+                basename_out = sww_file + '_' + quantity + '_' + extra_name_out
+
             file_out = sww2dem(dir+sep+sww_file, dir+sep+basename_out,
-                               quantity, 
+                               quantity,
                                timestep,
                                reduction,
@@ -1741 +1819 @@
                                format)
             files_out.append(file_out)
-    #print "basenames_out after",basenames_out 
     return files_out
 
 
+##
+# @brief 
+# @param production_dirs
+# @param output_dir
+# @param scenario_name
+# @param gauges_dir_name
+# @param plot_quantity
+# @param generate_fig
+# @param reportname
+# @param surface
+# @param time_min
+# @param time_max
+# @param title_on
+# @param verbose
+# @param nodes
 def get_timeseries(production_dirs, output_dir, scenario_name, gauges_dir_name,
-                   plot_quantity, generate_fig = False,
-                   reportname = None, surface = False, time_min = None,
-                   time_max = None, title_on = False, verbose = True,
+                   plot_quantity, generate_fig=False,
+                   reportname=None, surface=False, time_min=None,
+                   time_max=None, title_on=False, verbose=True,
                    nodes=None):
     """
@@ -1755 +1847 @@
     warning - this function has no tests
     """
+
     if reportname == None:
         report = False
     else:
         report = True
-        
+
     if nodes is None:
         is_parallel = False
     else:
         is_parallel = True
-        
+
     # Generate figures
     swwfiles = {}
-    
-    if is_parallel is True:    
+    if is_parallel is True:
         for i in range(nodes):
-            print 'Sending node %d of %d' %(i,nodes)
+            print 'Sending node %d of %d' % (i, nodes)
             swwfiles = {}
             if not reportname == None:
-                reportname = report_name + '_%s' %(i)
+                reportname = report_name + '_%s' % i
             for label_id in production_dirs.keys():
                 if label_id == 'boundaries':
@@ -1779 +1871 @@
                 else:
                     file_loc = output_dir + label_id + sep
-                    sww_extra = '_P%s_%s' %(i,nodes)
+                    sww_extra = '_P%s_%s' % (i, nodes)
                     swwfile = file_loc + scenario_name + sww_extra + '.sww'
-                    print 'swwfile',swwfile
+                    print 'swwfile', swwfile
                     swwfiles[swwfile] = label_id
 
@@ -1787 +1879 @@
                                               gauges_dir_name,
                                               production_dirs,
-                                              report = report,
-                                              reportname = reportname,
-                                              plot_quantity = plot_quantity,
-                                              generate_fig = generate_fig,
-                                              surface = surface,
-                                              time_min = time_min,
-                                              time_max = time_max,
-                                              title_on = title_on,
-                                              verbose = verbose)
-    else:   
-        for label_id in production_dirs.keys():       
+                                              report=report,
+                                              reportname=reportname,
+                                              plot_quantity=plot_quantity,
+                                              generate_fig=generate_fig,
+                                              surface=surface,
+                                              time_min=time_min,
+                                              time_max=time_max,
+                                              title_on=title_on,
+                                              verbose=verbose)
+    else:
+        for label_id in production_dirs.keys():
             if label_id == 'boundaries':
                 print 'boundaries'
@@ -1807 +1899 @@
                 swwfile = file_loc + scenario_name + '.sww'
             swwfiles[swwfile] = label_id
-        
+
         texname, elev_output = sww2timeseries(swwfiles,
                                               gauges_dir_name,
                                               production_dirs,
-                                              report = report,
-                                              reportname = reportname,
-                                              plot_quantity = plot_quantity,
-                                              generate_fig = generate_fig,
-                                              surface = surface,
-                                              time_min = time_min,
-                                              time_max = time_max,
-                                              title_on = title_on,
-                                              verbose = verbose)
-                                          
-
-    
-def sww2dem(basename_in, basename_out = None,
-            quantity = None, # defaults to elevation
-            timestep = None,
-            reduction = None,
-            cellsize = 10,
-            number_of_decimal_places = None,
-            NODATA_value = -9999,
-            easting_min = None,
-            easting_max = None,
-            northing_min = None,
-            northing_max = None,
-            verbose = False,
-            origin = None,
-            datum = 'WGS84',
-            format = 'ers'):
-
+                                              report=report,
+                                              reportname=reportname,
+                                              plot_quantity=plot_quantity,
+                                              generate_fig=generate_fig,
+                                              surface=surface,
+                                              time_min=time_min,
+                                              time_max=time_max,
+                                              title_on=title_on,
+                                              verbose=verbose)
+
+
+##
+# @brief Convert SWW file to DEM file (.asc or .ers).
+# @param basename_in 
+# @param basename_out 
+# @param quantity 
+# @param timestep 
+# @param reduction 
+# @param cellsize  
+# @param number_of_decimal_places 
+# @param NODATA_value 
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
+# @param verbose 
+# @param origin 
+# @param datum 
+# @param format 
+# @return 
+def sww2dem(basename_in, basename_out=None,
+            quantity=None, # defaults to elevation
+            timestep=None,
+            reduction=None,
+            cellsize=10,
+            number_of_decimal_places=None,
+            NODATA_value=-9999,
+            easting_min=None,
+            easting_max=None,
+            northing_min=None,
+            northing_max=None,
+            verbose=False,
+            origin=None,
+            datum='WGS84',
+            format='ers'):
     """Read SWW file and convert to Digitial Elevation model format
     (.asc or .ers)
 
     Example (ASC):
-
     ncols         3121
     nrows         1800
@@ -1854 +1962 @@
     The number of decimal places can be specified by the user to save
     on disk space requirements by specifying in the call to sww2dem.
-    
+
     Also write accompanying file with same basename_in but extension .prj
     used to fix the UTM zone, datum, false northings and eastings.
@@ -1885 +1993 @@
 
     import sys
-    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, \
-         sometrue
-    from Numeric import array2string
+    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape
+    from Numeric import array2string, sometrue
 
     from anuga.utilities.polygon import inside_polygon, outside_polygon, \
@@ -1897 +2004 @@
     assert format.lower() in ['asc', 'ers'], msg
 
-
     false_easting = 500000
     false_northing = 10000000
@@ -1903 +2009 @@
     if quantity is None:
         quantity = 'elevation'
-        
+
     if reduction is None:
         reduction = max
 
     if basename_out is None:
-        basename_out = basename_in + '_%s' %quantity
+        basename_out = basename_in + '_%s' % quantity
 
     if quantity_formula.has_key(quantity):
@@ -1915 +2021 @@
     if number_of_decimal_places is None:
         number_of_decimal_places = 3
-        
+
     swwfile = basename_in + '.sww'
     demfile = basename_out + '.' + format
     # Note the use of a .ers extension is optional (write_ermapper_grid will
     # deal with either option
-    
-    #if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors
-    
+
     # Read sww file
-    if verbose: 
-        print 'Reading from %s' %swwfile
-        print 'Output directory is %s' %basename_out
-    
+    if verbose:
+        print 'Reading from %s' % swwfile
+        print 'Output directory is %s' % basename_out
+
     from Scientific.IO.NetCDF import NetCDFFile
     fid = NetCDFFile(swwfile)
@@ -1942 +2046 @@
     number_of_timesteps = fid.dimensions['number_of_timesteps']
     number_of_points = fid.dimensions['number_of_points']
-    
+
     if origin is None:
-
         # Get geo_reference
         # sww files don't have to have a geo_ref
@@ -1959 +2062 @@
         xllcorner = origin[1]
         yllcorner = origin[2]
-
-
 
     # FIXME: Refactor using code from Interpolation_function.statistics
@@ -1996 +2097 @@
             q = fid.variables[name][:].flat
             if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
-            
+
     # Get quantity and reduce if applicable
     if verbose: print 'Processing quantity %s' %quantity
@@ -2002 +2103 @@
     # Turn NetCDF objects into Numeric arrays
     try:
-        q = fid.variables[quantity][:] 
-        
-        
+        q = fid.variables[quantity][:]
     except:
         quantity_dict = {}
         for name in fid.variables.keys():
-            quantity_dict[name] = fid.variables[name][:] 
-        #Convert quantity expression to quantities found in sww file    
+            quantity_dict[name] = fid.variables[name][:]
+        #Convert quantity expression to quantities found in sww file
         q = apply_expression_to_dictionary(quantity, quantity_dict)
-    #print "q.shape",q.shape
+
     if len(q.shape) == 2:
-        #q has a time component and needs to be reduced along
-        #the temporal dimension
+        #q has a time component, must be reduced alongthe temporal dimension
         if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = zeros( number_of_points, Float )
-        
+        q_reduced = zeros(number_of_points, Float)
+
         if timestep is not None:
             for k in range(number_of_points):
-                q_reduced[k] = q[timestep,k] 
+                q_reduced[k] = q[timestep,k]
         else:
             for k in range(number_of_points):
@@ -2032 +2130 @@
 
     if verbose:
-        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
-
+        print 'Processed values for %s are in [%f, %f]' % \
+              (quantity, min(q), max(q))
 
     #Create grid and update xll/yll corner and x,y
-
     #Relative extent
     if easting_min is None:
@@ -2057 +2154 @@
     else:
         ymax = northing_max - yllcorner
-    
 
     msg = 'xmax must be greater than or equal to xmin.\n'
@@ -2065 +2161 @@
     msg = 'yax must be greater than or equal to xmin.\n'
     msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
-    assert ymax >= ymin, msg    
-    
+    assert ymax >= ymin, msg
+
     if verbose: print 'Creating grid'
-    ncols = int((xmax-xmin)/cellsize)+1
-    nrows = int((ymax-ymin)/cellsize)+1
-
+    ncols = int((xmax-xmin)/cellsize) + 1
+    nrows = int((ymax-ymin)/cellsize) + 1
 
     #New absolute reference and coordinates
-    newxllcorner = xmin+xllcorner
-    newyllcorner = ymin+yllcorner
-
-    x = x+xllcorner-newxllcorner
-    y = y+yllcorner-newyllcorner
-    
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    newxllcorner = xmin + xllcorner
+    newyllcorner = ymin + yllcorner
+
+    x = x + xllcorner - newxllcorner
+    y = y + yllcorner - newyllcorner
+
+    vertex_points = concatenate ((x[:,NewAxis], y[:,NewAxis]), axis=1)
     assert len(vertex_points.shape) == 2
 
-    grid_points = zeros ( (ncols*nrows, 2), Float )
-
+    grid_points = zeros ((ncols*nrows, 2), Float)
 
     for i in xrange(nrows):
         if format.lower() == 'asc':
-            yg = i*cellsize
+            yg = i * cellsize
         else:
         #this will flip the order of the y values for ers
-            yg = (nrows-i)*cellsize
+            yg = (nrows-i) * cellsize
 
         for j in xrange(ncols):
-            xg = j*cellsize
+            xg = j * cellsize
             k = i*ncols + j
 
-            grid_points[k,0] = xg
-            grid_points[k,1] = yg
+            grid_points[k, 0] = xg
+            grid_points[k, 1] = yg
 
     #Interpolate
@@ -2112 +2206 @@
     grid_values = interp.interpolate(q, grid_points).flat
 
-
     if verbose:
         print 'Interpolated values are in [%f, %f]' %(min(grid_values),
@@ -2124 +2217 @@
         grid_values[i] = NODATA_value
 
-
-
-
     if format.lower() == 'ers':
         # setup ERS header information
-        grid_values = reshape(grid_values,(nrows, ncols))
+        grid_values = reshape(grid_values, (nrows, ncols))
         header = {}
         header['datum'] = '"' + datum + '"'
@@ -2148 +2238 @@
         #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
 
-
         #Write
         if verbose: print 'Writing %s' %demfile
+
         import ermapper_grids
+
         ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
 
@@ -2157 +2248 @@
     else:
         #Write to Ascii format
-
         #Write prj file
         prjfile = basename_out + '.prj'
@@ -2174 +2264 @@
         prjid.close()
 
-
-
         if verbose: print 'Writing %s' %demfile
 
@@ -2187 +2275 @@
         ascid.write('NODATA_value  %d\n' %NODATA_value)
 
-
         #Get bounding polygon from mesh
         #P = interp.mesh.get_boundary_polygon()
         #inside_indices = inside_polygon(grid_points, P)
-
         for i in range(nrows):
-            if verbose and i%((nrows+10)/10)==0:
+            if verbose and i % ((nrows+10)/10) == 0:
                 print 'Doing row %d of %d' %(i, nrows)
 
@@ -2200 +2286 @@
             slice = grid_values[base_index:base_index+ncols]
             #s = array2string(slice, max_line_width=sys.maxint)
-            s = array2string(slice, max_line_width=sys.maxint, precision=number_of_decimal_places)
+            s = array2string(slice, max_line_width=sys.maxint,
+                             precision=number_of_decimal_places)
             ascid.write(s[1:-1] + '\n')
-
-            #print
-            #for j in range(ncols):
-            #    index = base_index+j#
-            #    print grid_values[index],
-            #    ascid.write('%f ' %grid_values[index])
-            #ascid.write('\n')
 
         #Close
         ascid.close()
         fid.close()
+
         return basename_out
 
-
-#Backwards compatibility
+################################################################################
+# Obsolete functions - mainatined for backwards compatibility
+################################################################################
+
+##
+# @brief 
+# @param basename_in 
+# @param basename_out 
+# @param quantity 
+# @param timestep 
+# @param reduction 
+# @param cellsize 
+# @param verbose 
+# @param origin 
+# @note OBSOLETE - use sww2dem() instead.
 def sww2asc(basename_in, basename_out = None,
             quantity = None,
@@ -2237 +2331 @@
         format = 'asc')
 
-def sww2ers(basename_in, basename_out = None,
-            quantity = None,
-            timestep = None,
-            reduction = None,
-            cellsize = 10,
-            verbose = False,
-            origin = None,
-            datum = 'WGS84'):
+
+##
+# @brief 
+# @param basename_in
+# @param basename_out
+# @param quantity
+# @param timestep
+# @param reduction
+# @param cellsize
+# @param verbose
+# @param origin
+# @param datum
+# @note OBSOLETE - use sww2dem() instead.
+def sww2ers(basename_in, basename_out=None,
+            quantity=None,
+            timestep=None,
+            reduction=None,
+            cellsize=10,
+            verbose=False,
+            origin=None,
+            datum='WGS84'):
     print 'sww2ers will soon be obsoleted - please use sww2dem'
     sww2dem(basename_in,
-            basename_out = basename_out,
-            quantity = quantity,
-            timestep = timestep,
-            reduction = reduction,
-            cellsize = cellsize,
-            number_of_decimal_places = number_of_decimal_places,
-            verbose = verbose,
-            origin = origin,
-            datum = datum,
-            format = 'ers')
-################################# END COMPATIBILITY ##############
-
-
-
+            basename_out=basename_out,
+            quantity=quantity,
+            timestep=timestep,
+            reduction=reduction,
+            cellsize=cellsize,
+            number_of_decimal_places=number_of_decimal_places,
+            verbose=verbose,
+            origin=origin,
+            datum=datum,
+            format='ers')
+
+################################################################################
+# End of obsolete functions
+################################################################################
+
+
+##
+# @brief Convert SWW file to PTS file (at selected points).
+# @param basename_in Stem name of input SWW file.
+# @param basename_out Stem name of output file.
+# @param data_points If given, points where quantity is to be computed.
+# @param quantity Name (or expression) of existing quantity(s) (def: elevation).
+# @param timestep If given, output quantity at that timestep.
+# @param reduction If given, reduce quantity by this factor.
+# @param NODATA_value The NODATA value (default -9999).
+# @param verbose True if this function is to be verbose.
+# @param origin ??
 def sww2pts(basename_in, basename_out=None,
             data_points=None,
@@ -2269 +2389 @@
             verbose=False,
             origin=None):
-            #datum = 'WGS84')
-
-
     """Read SWW file and convert to interpolated values at selected points
 
-    The parameter quantity' must be the name of an existing quantity or
-    an expression involving existing quantities. The default is
-    'elevation'.
-
-    if timestep (an index) is given, output quantity at that timestep
+    The parameter 'quantity' must be the name of an existing quantity or
+    an expression involving existing quantities. The default is 'elevation'.
+
+    if timestep (an index) is given, output quantity at that timestep.
 
     if reduction is given use that to reduce quantity over all timesteps.
 
-    data_points (Nx2 array) give locations of points where quantity is to be computed.
-    
+    data_points (Nx2 array) give locations of points where quantity is to 
+    be computed.
     """
 
     import sys
-    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
-    from Numeric import array2string
-
-    from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
-    from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary
-
+    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape
+    from Numeric import array2string, sometrue
+    from anuga.utilities.polygon import inside_polygon, outside_polygon, \
+             separate_points_by_polygon
+    from anuga.abstract_2d_finite_volumes.util import \
+             apply_expression_to_dictionary
     from anuga.geospatial_data.geospatial_data import Geospatial_data
 
@@ -2302 +2418 @@
 
     if basename_out is None:
-        basename_out = basename_in + '_%s' %quantity
+        basename_out = basename_in + '_%s' % quantity
 
     swwfile = basename_in + '.sww'
@@ -2308 +2424 @@
 
     # Read sww file
-    if verbose: print 'Reading from %s' %swwfile
+    if verbose: print 'Reading from %s' % swwfile
     from Scientific.IO.NetCDF import NetCDFFile
     fid = NetCDFFile(swwfile)
@@ -2320 +2436 @@
     number_of_points = fid.dimensions['number_of_points']
     if origin is None:
-
         # Get geo_reference
         # sww files don't have to have a geo_ref
@@ -2326 +2441 @@
             geo_reference = Geo_reference(NetCDFObject=fid)
         except AttributeError, e:
-            geo_reference = Geo_reference() #Default georef object
+            geo_reference = Geo_reference() # Default georef object
 
         xllcorner = geo_reference.get_xllcorner()
@@ -2335 +2450 @@
         xllcorner = origin[1]
         yllcorner = origin[2]
-
-
 
     # FIXME: Refactor using code from file_function.statistics
@@ -2346 +2459 @@
         print '------------------------------------------------'
         print 'Statistics of SWW file:'
-        print '  Name: %s' %swwfile
+        print '  Name: %s' % swwfile
         print '  Reference:'
-        print '    Lower left corner: [%f, %f]'\
-              %(xllcorner, yllcorner)
-        print '    Start time: %f' %fid.starttime[0]
+        print '    Lower left corner: [%f, %f]' % (xllcorner, yllcorner)
+        print '    Start time: %f' % fid.starttime[0]
         print '  Extent:'
-        print '    x [m] in [%f, %f], len(x) == %d'\
-              %(min(x.flat), max(x.flat), len(x.flat))
-        print '    y [m] in [%f, %f], len(y) == %d'\
-              %(min(y.flat), max(y.flat), len(y.flat))
-        print '    t [s] in [%f, %f], len(t) == %d'\
-              %(min(times), max(times), len(times))
+        print '    x [m] in [%f, %f], len(x) == %d' \
+              % (min(x.flat), max(x.flat), len(x.flat))
+        print '    y [m] in [%f, %f], len(y) == %d' \
+              % (min(y.flat), max(y.flat), len(y.flat))
+        print '    t [s] in [%f, %f], len(t) == %d' \
+              % (min(times), max(times), len(times))
         print '  Quantities [SI units]:'
         for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
             q = fid.variables[name][:].flat
-            print '    %s in [%f, %f]' %(name, min(q), max(q))
-
-
+            print '    %s in [%f, %f]' % (name, min(q), max(q))
 
     # Get quantity and reduce if applicable
-    if verbose: print 'Processing quantity %s' %quantity
+    if verbose: print 'Processing quantity %s' % quantity
 
     # Turn NetCDF objects into Numeric arrays
@@ -2373 +2483 @@
         quantity_dict[name] = fid.variables[name][:]
 
-
-
-    # Convert quantity expression to quantities found in sww file    
+    # Convert quantity expression to quantities found in sww file
     q = apply_expression_to_dictionary(quantity, quantity_dict)
-
-
 
     if len(q.shape) == 2:
         # q has a time component and needs to be reduced along
         # the temporal dimension
-        if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = zeros( number_of_points, Float )
-
+        if verbose: print 'Reducing quantity %s' % quantity
+
+        q_reduced = zeros(number_of_points, Float)
         for k in range(number_of_points):
-            q_reduced[k] = reduction( q[:,k] )
-
+            q_reduced[k] = reduction(q[:,k])
         q = q_reduced
 
@@ -2395 +2500 @@
     assert q.shape[0] == number_of_points
 
-
     if verbose:
-        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
-
+        print 'Processed values for %s are in [%f, %f]' \
+              % (quantity, min(q), max(q))
 
     # Create grid and update xll/yll corner and x,y
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    vertex_points = concatenate((x[:, NewAxis], y[:, NewAxis]), axis=1)
     assert len(vertex_points.shape) == 2
 
     # Interpolate
     from anuga.fit_interpolate.interpolate import Interpolate
-    interp = Interpolate(vertex_points, volumes, verbose = verbose)
+    interp = Interpolate(vertex_points, volumes, verbose=verbose)
 
     # Interpolate using quantity values
@@ -2412 +2516 @@
     interpolated_values = interp.interpolate(q, data_points).flat
 
-
     if verbose:
-        print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
-                                                      max(interpolated_values))
-
-    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
+        print 'Interpolated values are in [%f, %f]' % (min(interpolated_values),
+                                                       max(interpolated_values))
+
+    # Assign NODATA_value to all points outside bounding polygon
+    # (from interpolation mesh)
     P = interp.mesh.get_boundary_polygon()
     outside_indices = outside_polygon(data_points, P, closed=True)
@@ -2424 +2528 @@
         interpolated_values[i] = NODATA_value
 
-
-    # Store results    
-    G = Geospatial_data(data_points=data_points,
-                        attributes=interpolated_values)
+    # Store results
+    G = Geospatial_data(data_points=data_points, attributes=interpolated_values)
 
     G.export_points_file(ptsfile, absolute = True)
@@ -2434 +2536 @@
 
 
-def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
-                                  use_cache = False,
-                                  verbose = False):
-    """Read Digitial Elevation model from the following ASCII format (.asc)
+##
+# @brief Convert ASC file to DEM file.
+# @param basename_in Stem of input filename.
+# @param basename_out Stem of output filename.
+# @param use_cache ??
+# @param verbose True if this function is to be verbose.
+# @return 
+# @note A PRJ file with same stem basename must exist and is used to fix the
+#       UTM zone, datum, false northings and eastings.
+def convert_dem_from_ascii2netcdf(basename_in, basename_out=None,
+                                  use_cache=False,
+                                  verbose=False):
+    """Read Digital Elevation model from the following ASCII format (.asc)
 
     Example:
-
     ncols         3121
     nrows         1800
@@ -2451 +2561 @@
     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
@@ -2469 +2578 @@
     """
 
-
-
     kwargs = {'basename_out': basename_out, 'verbose': verbose}
 
@@ -2476 +2583 @@
         from caching import cache
         result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
-                       dependencies = [basename_in + '.asc',
-                                       basename_in + '.prj'],
-                       verbose = verbose)
+                       dependencies=[basename_in + '.asc',
+                                     basename_in + '.prj'],
+                       verbose=verbose)
 
     else:
@@ -2486 +2593 @@
 
 
-
-
-
-
+##
+# @brief Convert an ASC file to a DEM file.
+# @param basename_in Stem of input filename.
+# @param basename_out Stem of output filename.
+# @param verbose True if this function is to be verbose.
 def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
                                   verbose = False):
-    """Read Digitial Elevation model from the following ASCII format (.asc)
-
-    Internal function. See public function convert_dem_from_ascii2netcdf for details.
+    """Read Digital Elevation model from the following ASCII format (.asc)
+
+    Internal function. See public function convert_dem_from_ascii2netcdf
+    for details.
     """
 
@@ -2501 +2610 @@
     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'
+    if verbose: print 'Reading METADATA from %s' % root + '.prj'
+
     metadatafile = open(root + '.prj')
     metalines = metadatafile.readlines()
@@ -2543 +2651 @@
     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')
+    if verbose: print 'Reading DEM from %s' % basename_in + '.asc'
+
     lines = datafile.readlines()
     datafile.close()
@@ -2561 +2665 @@
 
     # Do cellsize (line 4) before line 2 and 3
-    cellsize = float(lines[4].split()[1].strip())    
+    cellsize = float(lines[4].split()[1].strip())
 
     # Checks suggested by Joaquim Luis
@@ -2572 +2676 @@
         xllcorner = float(xref[1].strip())
     else:
-        msg = 'Unknown keyword: %s' %xref[0].strip()
+        msg = 'Unknown keyword: %s' % xref[0].strip()
         raise Exception, msg
 
@@ -2581 +2685 @@
         yllcorner = float(yref[1].strip())
     else:
-        msg = 'Unknown keyword: %s' %yref[0].strip()
+        msg = 'Unknown keyword: %s' % yref[0].strip()
         raise Exception, msg
-        
 
     NODATA_value = int(lines[5].split()[1].strip())
 
     assert len(lines) == nrows + 6
-
-
-    ##########################################
-
 
     if basename_out == None:
@@ -2598 +2697 @@
         netcdfname = basename_out + '.dem'
 
-    if verbose: print 'Store to NetCDF file %s' %netcdfname
+    if verbose: print 'Store to NetCDF file %s' % netcdfname
+
     # NetCDF file definition
     fid = NetCDFFile(netcdfname, 'w')
@@ -2604 +2704 @@
     #Create new file
     fid.institution = 'Geoscience Australia'
-    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
+    fid.description = 'NetCDF DEM format for compact and portable storage ' \
                       'of spatial point data'
 
@@ -2621 +2721 @@
     fid.units = units
 
-
     # dimension definitions
     fid.createDimension('number_of_rows', nrows)
@@ -2637 +2736 @@
     for i, line in enumerate(lines[6:]):
         fields = line.split()
-        if verbose and i%((n+10)/10)==0:
-            print 'Processing row %d of %d' %(i, nrows)
-
+        if verbose and i % ((n+10)/10) == 0:
+            print 'Processing row %d of %d' % (i, nrows)
         elevation[i, :] = array([float(x) for x in fields])
 
@@ -2645 +2743 @@
 
 
-
-
-
-def ferret2sww(basename_in, basename_out = None,
-               verbose = False,
-               minlat = None, maxlat = None,
-               minlon = None, maxlon = None,
-               mint = None, maxt = None, mean_stage = 0,
-               origin = None, zscale = 1,
-               fail_on_NaN = True,
-               NaN_filler = 0,
-               elevation = None,
-               inverted_bathymetry = True
+##
+# @brief Convert 'ferret' file to SWW file.
+# @param basename_in Stem of input filename.
+# @param basename_out Stem of output filename.
+# @param verbose True if this function is to be verbose.
+# @param minlat 
+# @param maxlat 
+# @param minlon 
+# @param maxlon 
+# @param mint 
+# @param maxt 
+# @param mean_stage 
+# @param origin 
+# @param zscale 
+# @param fail_on_NaN 
+# @param NaN_filler 
+# @param elevation 
+# @param inverted_bathymetry 
+def ferret2sww(basename_in, basename_out=None,
+               verbose=False,
+               minlat=None, maxlat=None,
+               minlon=None, maxlon=None,
+               mint=None, maxt=None, mean_stage=0,
+               origin=None, zscale=1,
+               fail_on_NaN=True,
+               NaN_filler=0,
+               elevation=None,
+               inverted_bathymetry=True
                ): #FIXME: Bathymetry should be obtained
                                   #from MOST somehow.
@@ -2712 +2825 @@
         if minlon != None:
             assert maxlon > minlon
-        
-
-
-    #Get NetCDF data
-    if verbose: print 'Reading files %s_*.nc' %basename_in
-    #print "basename_in + '_ha.nc'",basename_in + '_ha.nc' 
-    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
-    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
-    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
-    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
+
+    # Get NetCDF data
+    if verbose: print 'Reading files %s_*.nc' % basename_in
+
+    file_h = NetCDFFile(basename_in + '_ha.nc', 'r') # Wave amplitude (cm)
+    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') # Velocity (x) (cm/s)
+    file_v = NetCDFFile(basename_in + '_va.nc', 'r') # Velocity (y) (cm/s)
+    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  # Elevation (z) (m)
 
     if basename_out is None:
@@ -2736 +2847 @@
         if dimension[:4] == 'TIME':
             dim_h_time = dimension
-
-#    print 'long:', dim_h_longitude
-#    print 'lats:', dim_h_latitude
-#    print 'times:', dim_h_time
 
     times = file_h.variables[dim_h_time]
@@ -2764 +2871 @@
         if dimension[:4] == 'TIME':
             dim_u_time = dimension
-            
+
     # get dimensions for file_v
     for dimension in file_v.dimensions.keys():
@@ -2774 +2881 @@
             dim_v_time = dimension
 
-
     # Precision used by most for lat/lon is 4 or 5 decimals
     e_lat = around(file_e.variables[dim_e_latitude][:], 5)
@@ -2790 +2896 @@
     if mint is None:
         jmin = 0
-        mint = times[0]        
+        mint = times[0]
     else:
         jmin = searchsorted(times, mint)
-        
+
     if maxt is None:
         jmax = len(times)
@@ -2800 +2906 @@
         jmax = searchsorted(times, maxt)
 
-    #print "latitudes[:]",latitudes[:]
-    #print "longitudes[:]",longitudes [:]
     kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
                                                   longitudes[:],
@@ -2812 +2916 @@
     longitudes = longitudes[lmin:lmax]
 
-    #print "latitudes[:]",latitudes[:]
-    #print "longitudes[:]",longitudes [:]
-
     if verbose: print 'cropping'
+
     zname = 'ELEVATION'
 
@@ -2837 +2939 @@
     #        'print hmmm'
 
-
-
     #Get missing values
     nan_ha = file_h.variables['HA'].missing_value[0]
@@ -2854 +2954 @@
     if sometrue (missing):
         if fail_on_NaN:
-            msg = 'NetCDFFile %s contains missing values'\
-                  %(basename_in+'_ha.nc')
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_ha.nc'
             raise DataMissingValuesError, msg
         else:
@@ -2863 +2963 @@
     if sometrue (missing):
         if fail_on_NaN:
-            msg = 'NetCDFFile %s contains missing values'\
-                  %(basename_in+'_ua.nc')
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_ua.nc'
             raise DataMissingValuesError, msg
         else:
@@ -2872 +2972 @@
     if sometrue (missing):
         if fail_on_NaN:
-            msg = 'NetCDFFile %s contains missing values'\
-                  %(basename_in+'_va.nc')
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_va.nc'
             raise DataMissingValuesError, msg
         else:
             vspeed = vspeed*(missing==0) + missing*NaN_filler
 
-
     missing = (elevations == nan_e)
     if sometrue (missing):
         if fail_on_NaN:
-            msg = 'NetCDFFile %s contains missing values'\
-                  %(basename_in+'_e.nc')
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_e.nc'
             raise DataMissingValuesError, msg
         else:
             elevations = elevations*(missing==0) + missing*NaN_filler
-
-    #######
-
-
 
     number_of_times = times.shape[0]
@@ -2904 +2999 @@
         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))
+        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))
+        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))
+        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 '  %s in [%f, %f]' % (name, min(q), max(q))
 
         q = elevations.flat
         name = 'Elevations (e) [m]'
-        print '  %s in [%f, %f]' %(name, min(q), max(q))
-
+        print '  %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
-
+    number_of_points = number_of_latitudes * number_of_longitudes
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
 
     file_h.close()
@@ -2940 +3032 @@
     file_e.close()
 
-
     # NetCDF file definition
     outfile = NetCDFFile(swwname, 'w')
 
-    description = 'Converted from Ferret files: %s, %s, %s, %s'\
-                  %(basename_in + '_ha.nc',
-                    basename_in + '_ua.nc',
-                    basename_in + '_va.nc',
-                    basename_in + '_e.nc')
-    
+    description = 'Converted from Ferret files: %s, %s, %s, %s' \
+                  % (basename_in + '_ha.nc',
+                     basename_in + '_ua.nc',
+                     basename_in + '_va.nc',
+                     basename_in + '_e.nc')
+
     # Create new file
     starttime = times[0]
-    
+
     sww = Write_sww()
     sww.store_header(outfile, times, number_of_volumes,
                      number_of_points, description=description,
-                     verbose=verbose,sww_precision=Float)
+                     verbose=verbose, sww_precision=Float)
 
     # Store
@@ -2963 +3054 @@
     y = zeros(number_of_points, Float)  #Northing
 
-
     if verbose: print 'Making triangular grid'
 
     # Check zone boundaries
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+    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
-
+    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
+            #msg = 'Zone boundary crossed at longitude =', lon
             #assert zone == refzone, msg
             #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
@@ -2987 +3076 @@
     #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
+    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]
@@ -3000 +3089 @@
 
     if origin is None:
-        origin = Geo_reference(refzone,min(x),min(y))
+        origin = Geo_reference(refzone, min(x), min(y))
     geo_ref = write_NetCDF_georeference(origin, outfile)
-    
+
     if elevation is not None:
         z = elevation
     else:
         if inverted_bathymetry:
-            z = -1*elevations
+            z = -1 * elevations
         else:
             z = elevations
@@ -3014 +3103 @@
     #FIXME use the Write_sww instance(sww) to write this info
     from Numeric import resize
-    z = resize(z,outfile.variables['z'][:].shape)
+    z = resize(z, outfile.variables['z'][:].shape)
     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
@@ -3021 +3110 @@
     outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
 
-
-
     #Time stepping
     stage = outfile.variables['stage']
@@ -3029 +3116 @@
 
     if verbose: print 'Converting quantities'
+
     n = len(times)
     for j in range(n):
-        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
+        if verbose and j % ((n+10)/10) == 0: print '  Doing %d of %d' %(j, n)
+
         i = 0
-        for k in range(number_of_latitudes):      #Y direction
-            for l in range(number_of_longitudes): #X direction
-                w = zscale*amplitudes[j,k,l]/100 + mean_stage
+        for k in range(number_of_latitudes):      # Y direction
+            for l in range(number_of_longitudes): # X direction
+                w = zscale * amplitudes[j,k,l] / 100 + mean_stage
                 stage[j,i] = w
                 h = w - z[i]
@@ -3053 +3142 @@
         print '  Name: %s' %swwname
         print '  Reference:'
-        print '    Lower left corner: [%f, %f]'\
-              %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
+        print '    Lower left corner: [%f, %f]' \
+              % (geo_ref.get_xllcorner(), geo_ref.get_yllcorner())
         print ' Start time: %f' %starttime
         print '    Min time: %f' %mint
         print '    Max time: %f' %maxt
         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 '    x [m] in [%f, %f], len(x) == %d' \
+              % (min(x.flat), max(x.flat), len(x.flat))
+        print '    y [m] in [%f, %f], len(y) == %d' \
+              % (min(y.flat), max(y.flat), len(y.flat))
+        print '    t [s] in [%f, %f], len(t) == %d' \
+              % (min(times), max(times), len(times))
         print '  Quantities [SI units]:'
         for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
             q = outfile.variables[name][:].flat
-            print '    %s in [%f, %f]' %(name, min(q), max(q))
-
-
+            print '    %s in [%f, %f]' % (name, min(q), max(q))
 
     outfile.close()
 
 
-
-
-
+##
+# @brief Convert time-series text file to TMS file.
+# @param filename 
+# @param quantity_names 
+# @param time_as_seconds 
 def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
     """Template for converting typical text files with time series to
     NetCDF tms file.
 
-
     The file format is assumed to be either two fields separated by a comma:
 
@@ -3093 +3181 @@
 
     or time (seconds), value0 value1 value2 ...
-    
+
       0.0, 1.328223 0 0
       0.1, 1.292912 0 0
@@ -3107 +3195 @@
     from anuga.utilities.numerical_tools import ensure_numeric
 
-
-
-    fid = open(filename + '.txt')
+    file_text = filename + '.txt'
+    fid = open(file_text)
     line = fid.readline()
     fid.close()
 
     fields = line.split(',')
-    msg = 'File %s must have the format date, value0 value1 value2 ...'
+    msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \
+          % file_text
     assert len(fields) == 2, msg
 
@@ -3121 +3209 @@
             starttime = calendar.timegm(time.strptime(fields[0], time_format))
         except ValueError:
-            msg = 'First field in file %s must be' %filename
-            msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0]
+            msg = 'First field in file %s must be' % file_text
+            msg += ' date-time with format %s.\n' % time_format
+            msg += 'I got %s instead.' % fields[0]
             raise DataTimeError, msg
     else:
@@ -3132 +3220 @@
             raise DataTimeError, msg
 
-
-    #Split values
+    # Split values
     values = []
     for value in fields[1].split():
@@ -3143 +3230 @@
     assert len(q.shape) == 1, msg
 
-
-
-    #Read times proper
+    # Read times proper
     from Numeric import zeros, Float, alltrue
     from anuga.config import time_format
     import time, calendar
 
-    fid = open(filename + '.txt')
+    fid = open(file_text)
     lines = fid.readlines()
     fid.close()
@@ -3157 +3242 @@
     d = len(q)
 
-    T = zeros(N, Float)       #Time
-    Q = zeros((N, d), Float)  #Values
+    T = zeros(N, Float)       # Time
+    Q = zeros((N, d), Float)  # Values
 
     for i, line in enumerate(lines):
@@ -3171 +3256 @@
             Q[i, j] = float(value)
 
-    msg = 'File %s must list time as a monotonuosly ' %filename
+    msg = 'File %s must list time as a monotonuosly ' % filename
     msg += 'increasing sequence'
-    assert alltrue( T[1:] - T[:-1] > 0 ), msg
+    assert alltrue(T[1:] - T[:-1] > 0), msg
 
     #Create NetCDF file
@@ -3180 +3265 @@
     fid = NetCDFFile(filename + '.tms', 'w')
 
-
     fid.institution = 'Geoscience Australia'
     fid.description = 'Time series'
-
 
     #Reference point
@@ -3194 +3277 @@
     #fid.createDimension('number_of_vertices', 3)
 
-
     fid.createDimension('number_of_timesteps', len(T))
 
@@ -3205 +3287 @@
             name = quantity_names[i]
         except:
-            name = 'Attribute%d'%i
+            name = 'Attribute%d' % i
 
         fid.createVariable(name, Float, ('number_of_timesteps',))
@@ -3213 +3295 @@
 
 
+##
+# @brief Get the extents of a NetCDF data file.
+# @param file_name The path to the SWW file.
+# @return A list of x, y, z and stage limits (min, max).
 def extent_sww(file_name):
-    """
-    Read in an sww file.
-
-    Input;
+    """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)
-    """
-
+    Output:
+    A list: [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
+    """
 
     from Scientific.IO.NetCDF import NetCDFFile
 
-    #Check contents
     #Get NetCDF
     fid = NetCDFFile(file_name, 'r')
@@ -3239 +3318 @@
     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)]
-
-
+    return [min(x), max(x), min(y), max(y), min(stage.flat), max(stage.flat)]
+
+
+##
+# @brief 
+# @param filename
+# @param boundary
+# @param t
+# @param fail_if_NaN
+# @param NaN_filler
+# @param verbose
+# @param very_verbose
+# @return 
 def sww2domain(filename, boundary=None, t=None,
-               fail_if_NaN=True ,NaN_filler=0,
-               verbose = False, very_verbose = False):
+               fail_if_NaN=True, NaN_filler=0,
+               verbose=False, very_verbose=False):
     """
     Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
@@ -3260 +3345 @@
     give a different final boundary, or crash.
     """
-    NaN=9.969209968386869e+036
-    #initialise NaN.
 
     from Scientific.IO.NetCDF import NetCDFFile
@@ -3267 +3350 @@
     from Numeric import asarray, transpose, resize
 
+    # initialise NaN.
+    NaN = 9.969209968386869e+036
+
     if verbose: print 'Reading from ', filename
-    fid = NetCDFFile(filename, 'r')    #Open existing file for read
-    time = fid.variables['time']       #Timesteps
+
+    fid = NetCDFFile(filename, 'r')    # Open existing file for read
+    time = fid.variables['time']       # Timesteps
     if t is None:
         t = time[-1]
@@ -3275 +3362 @@
 
     # 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
+    x = fid.variables['x'][:]                   # x-coordinates of vertices
+    y = fid.variables['y'][:]                   # y-coordinates of vertices
+    elevation = fid.variables['elevation']      # Elevation
+    stage = fid.variables['stage']              # Water level
+    xmomentum = fid.variables['xmomentum']      # Momentum in the x-direction
+    ymomentum = fid.variables['ymomentum']      # Momentum in the y-direction
 
     starttime = fid.starttime[0]
-    volumes = fid.variables['volumes'][:] #Connectivity
-    coordinates = transpose(asarray([x.tolist(),y.tolist()]))
-    #FIXME (Ole): Something like this might be better: concatenate( (x, y), axis=1 )
-    # or concatenate( (x[:,NewAxis],x[:,NewAxis]), axis=1 )
+    volumes = fid.variables['volumes'][:]       # Connectivity
+    coordinates = transpose(asarray([x.tolist(), y.tolist()]))
+    # FIXME (Ole): Something like this might be better:
+    #                 concatenate((x, y), axis=1)
+    # or              concatenate((x[:,NewAxis], x[:,NewAxis]), axis=1)
 
     conserved_quantities = []
@@ -3293 +3381 @@
 
     # get geo_reference
-    #sww files don't have to have a geo_ref
-    try:
+    try:                             # sww files don't have to have a geo_ref
         geo_reference = Geo_reference(NetCDFObject=fid)
-    except: #AttributeError, e:
+    except: # AttributeError, e:
         geo_reference = None
 
     if verbose: print '    getting quantities'
+
     for quantity in fid.variables.keys():
         dimensions = fid.variables[quantity].dimensions
         if 'number_of_timesteps' in dimensions:
             conserved_quantities.append(quantity)
-            interpolated_quantities[quantity]=\
-                  interpolated_quantity(fid.variables[quantity][:],time_interp)
-        else: other_quantities.append(quantity)
+            interpolated_quantities[quantity] = \
+                  interpolated_quantity(fid.variables[quantity][:], time_interp)
+        else:
+            other_quantities.append(quantity)
 
     other_quantities.remove('x')
@@ -3319 +3408 @@
     except:
         pass
-        
 
     conserved_quantities.remove('time')
 
     if verbose: print '    building domain'
+
     #    From domain.Domain:
     #    domain = Domain(coordinates, volumes,\
@@ -3330 +3419 @@
     #                    xllcorner=xllcorner, yllcorner=yllcorner)
 
-    #   From shallow_water.Domain:
-    coordinates=coordinates.tolist()
-    volumes=volumes.tolist()
-    #FIXME:should this be in mesh?(peter row)
-    if fid.smoothing == 'Yes': unique = False
-    else: unique = True
+    # From shallow_water.Domain:
+    coordinates = coordinates.tolist()
+    volumes = volumes.tolist()
+    # FIXME:should this be in mesh? (peter row)
+    if fid.smoothing == 'Yes':
+        unique = False
+    else:
+        unique = True
     if unique:
-        coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
-
+        coordinates, volumes, boundary = weed(coordinates, volumes,boundary)
 
     try:
@@ -3344 +3434 @@
     except AssertionError, e:
         fid.close()
-        msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e
+        msg = 'Domain could not be created: %s. ' \
+              'Perhaps use "fail_if_NaN=False and NaN_filler = ..."' % e
         raise DataDomainError, msg
 
@@ -3352 +3443 @@
     domain.geo_reference = geo_reference
 
-    domain.starttime=float(starttime)+float(t)
-    domain.time=0.0
+    domain.starttime = float(starttime) + float(t)
+    domain.time = 0.0
 
     for quantity in other_quantities:
@@ -3359 +3450 @@
             NaN = fid.variables[quantity].missing_value
         except:
-            pass #quantity has no missing_value number
+            pass                       # quantity has no missing_value number
         X = fid.variables[quantity][:]
         if very_verbose:
-            print '       ',quantity
-            print '        NaN =',NaN
+            print '       ', quantity
+            print '        NaN =', NaN
             print '        max(X)'
-            print '       ',max(X)
+            print '       ', max(X)
             print '        max(X)==NaN'
-            print '       ',max(X)==NaN
+            print '       ', max(X)==NaN
             print ''
-        if (max(X)==NaN) or (min(X)==NaN):
+        if max(X) == NaN or min(X) == NaN:
             if fail_if_NaN:
-                msg = 'quantity "%s" contains no_data entry'%quantity
+                msg = 'quantity "%s" contains no_data entry' % quantity
                 raise DataMissingValuesError, msg
             else:
-                data = (X<>NaN)
-                X = (X*data)+(data==0)*NaN_filler
+                data = (X != NaN)
+                X = (X*data) + (data==0)*NaN_filler
         if unique:
-            X = resize(X,(len(X)/3,3))
-        domain.set_quantity(quantity,X)
+            X = resize(X, (len(X)/3, 3))
+        domain.set_quantity(quantity, X)
     #
     for quantity in conserved_quantities:
@@ -3384 +3475 @@
             NaN = fid.variables[quantity].missing_value
         except:
-            pass #quantity has no missing_value number
+            pass                       # quantity has no missing_value number
         X = interpolated_quantities[quantity]
         if very_verbose:
             print '       ',quantity
-            print '        NaN =',NaN
+            print '        NaN =', NaN
             print '        max(X)'
-            print '       ',max(X)
+            print '       ', max(X)
             print '        max(X)==NaN'
-            print '       ',max(X)==NaN
+            print '       ', max(X)==NaN
             print ''
-        if (max(X)==NaN) or (min(X)==NaN):
+        if max(X) == NaN or min(X) == NaN:
             if fail_if_NaN:
-                msg = 'quantity "%s" contains no_data entry'%quantity
+                msg = 'quantity "%s" contains no_data entry' % quantity
                 raise DataMissingValuesError, msg
             else:
-                data = (X<>NaN)
-                X = (X*data)+(data==0)*NaN_filler
+                data = (X != NaN)
+                X = (X*data) + (data==0)*NaN_filler
         if unique:
-            X = resize(X,(X.shape[0]/3,3))
-        domain.set_quantity(quantity,X)
+            X = resize(X, (X.shape[0]/3, 3))
+        domain.set_quantity(quantity, X)
 
     fid.close()
+
     return domain
 
 
-def interpolated_quantity(saved_quantity,time_interp):
-
-    #given an index and ratio, interpolate quantity with respect to time.
-    index,ratio = time_interp
+##
+# @brief Interpolate a quantity wrt time.
+# @param saved_quantity The quantity to interpolate.
+# @param time_interp (index, ratio)
+# @return A vector of interpolated values.
+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]
+        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):
+##
+# @brief 
+# @parm time 
+# @param t 
+# @return An (index, ration) tuple.
+def get_time_interp(time, t=None):
     #Finds the ratio and index for time interpolation.
     #It is borrowed from previous abstract_2d_finite_volumes code.
@@ -3433 +3538 @@
         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
+        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 DataTimeError, msg
         if tau > time[-1]: raise DataTimeError, msg
@@ -3447 +3552 @@
             #t is now between index and index+1
             ratio = (tau - time[index])/(time[index+1] - time[index])
-    return (index,ratio)
-
-
-def weed(coordinates,volumes,boundary = None):
-    if type(coordinates)==ArrayType:
+
+    return (index, ratio)
+
+
+##
+# @brief 
+# @param coordinates 
+# @param volumes 
+# @param boundary 
+def weed(coordinates, volumes, boundary=None):
+    if type(coordinates) == ArrayType:
         coordinates = coordinates.tolist()
-    if type(volumes)==ArrayType:
+    if type(volumes) == ArrayType:
         volumes = volumes.tolist()
 
@@ -3463 +3574 @@
         if point_dict.has_key(point):
             unique = True
-            same_point[i]=point
+            same_point[i] = point
             #to change all point i references to point j
         else:
-            point_dict[point]=i
-            same_point[i]=point
+            point_dict[point] = i
+            same_point[i] = point
 
     coordinates = []
@@ -3474 +3585 @@
         point = tuple(point)
         coordinates.append(list(point))
-        point_dict[point]=i
-        i+=1
-
+        point_dict[point] = i
+        i += 1
 
     for volume in volumes:
         for i in range(len(volume)):
             index = volume[i]
-            if index>-1:
-                volume[i]=point_dict[same_point[index]]
+            if index > -1:
+                volume[i] = point_dict[same_point[index]]
 
     new_boundary = {}
@@ -3493 +3603 @@
             #('exterior, pond') or replaced ('pond')(peter row)
 
-            if new_boundary.has_key((point0,point1)):
-                new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\
-                                              #+','+label
-
-            elif new_boundary.has_key((point1,point0)):
-                new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\
-                                              #+','+label
-            else: new_boundary[(point0,point1)]=label
+            if new_boundary.has_key((point0, point1)):
+                new_boundary[(point0,point1)] = new_boundary[(point0, point1)]
+
+            elif new_boundary.has_key((point1, point0)):
+                new_boundary[(point1,point0)] = new_boundary[(point1, point0)]
+            else: new_boundary[(point0, point1)] = label
 
         boundary = new_boundary
 
-    return coordinates,volumes,boundary
-
-
+    return coordinates, volumes, boundary
+
+
+##
+# @brief Read DEM file, decimate data, write new DEM file.
+# @param basename_in Stem of the input filename.
+# @param stencil 
+# @param cellsize_new New cell size to resample on.
+# @param basename_out Stem of the output filename.
+# @param verbose True if this function is to be verbose.
 def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
                  verbose=False):
@@ -3533 +3648 @@
     #Open existing netcdf file to read
     infile = NetCDFFile(inname, 'r')
-    if verbose: print 'Reading DEM from %s' %inname
+
+    if verbose: print 'Reading DEM from %s' % inname
 
     #Read metadata
@@ -3557 +3673 @@
         outname = basename_out + '.dem'
 
-    if verbose: print 'Write decimated NetCDF file to %s' %outname
+    if verbose: print 'Write decimated NetCDF file to %s' % outname
 
     #Determine some dimensions for decimated grid
@@ -3572 +3688 @@
     #Create new file
     outfile.institution = 'Geoscience Australia'
-    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
-                           'of spatial point data'
+    outfile.description = 'NetCDF DEM format for compact and portable ' \
+                          'storage of spatial point data'
+
     #Georeferencing
     outfile.zone = zone
@@ -3605 +3722 @@
     for i in range(nrows_new):
         if verbose: print 'Processing row %d of %d' %(i, nrows_new)
+
         lower_index = global_index
-        telev =  zeros(ncols_new, Float)
+        telev = zeros(ncols_new, Float)
         local_index = 0
         trow = i * cellsize_ratio
@@ -3612 +3730 @@
         for j in range(ncols_new):
             tcol = j * cellsize_ratio
-            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
+            tmp = dem_elevation_r[trow:trow+nrows_stencil,
+                                  tcol:tcol+ncols_stencil]
 
             #if dem contains 1 or more NODATA_values set value in
@@ -3629 +3748 @@
         elevation[lower_index:upper_index] = telev
 
-    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
+    assert global_index == nrows_new*ncols_new, \
+           'index not equal to number of points'
 
     infile.close()
@@ -3635 +3755 @@
 
 
-
-
-def tsh2sww(filename, verbose=False): 
+##
+# @brief 
+# @param filename
+# @param verbose
+def tsh2sww(filename, verbose=False):
     """
     to check if a tsh/msh file 'looks' good.
     """
 
-
     if verbose == True:print 'Creating domain from', filename
+
     domain = pmesh_to_domain_instance(filename, Domain)
+
     if verbose == True:print "Number of triangles = ", len(domain)
 
@@ -3651 +3774 @@
     file_path, filename = path.split(filename)
     filename, ext = path.splitext(filename)
-    domain.set_name(filename)    
+    domain.set_name(filename)
     domain.reduction = mean
+
     if verbose == True:print "file_path",file_path
-    if file_path == "":file_path = "."
+
+    if file_path == "":
+        file_path = "."
     domain.set_datadir(file_path)
 
@@ -3660 +3786 @@
         print "Output written to " + domain.get_datadir() + sep + \
               domain.get_name() + "." + domain.format
+
     sww = get_dataobject(domain)
     sww.store_connectivity()
@@ -3665 +3792 @@
 
 
+##
+# @brief Convert CSIRO ESRI file to an SWW boundary file.
+# @param bath_dir 
+# @param elevation_dir 
+# @param ucur_dir 
+# @param vcur_dir 
+# @param sww_file 
+# @param minlat 
+# @param maxlat 
+# @param minlon 
+# @param maxlon 
+# @param zscale 
+# @param mean_stage 
+# @param fail_on_NaN 
+# @param elevation_NaN_filler 
+# @param bath_prefix 
+# @param elevation_prefix 
+# @param verbose 
+# @note Also convert latitude and longitude to UTM. All coordinates are
+#       assumed to be given in the GDA94 datum.
 def asc_csiro2sww(bath_dir,
                   elevation_dir,
@@ -3670 +3817 @@
                   vcur_dir,
                   sww_file,
-                  minlat = None, maxlat = None,
-                  minlon = None, maxlon = None,
+                  minlat=None, maxlat=None,
+                  minlon=None, maxlon=None,
                   zscale=1,
-                  mean_stage = 0,
-                  fail_on_NaN = True,
-                  elevation_NaN_filler = 0,
+                  mean_stage=0,
+                  fail_on_NaN=True,
+                  elevation_NaN_filler=0,
                   bath_prefix='ba',
                   elevation_prefix='el',
@@ -3700 +3847 @@
     The time period is less than 24hrs and uniform.
     """
+
     from Scientific.IO.NetCDF import NetCDFFile
 
@@ -3708 +3856 @@
     # go in to the bath dir and load the only file,
     bath_files = os.listdir(bath_dir)
-
     bath_file = bath_files[0]
     bath_dir_file =  bath_dir + os.sep + bath_file
-    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
+    bath_metadata, bath_grid =  _read_asc(bath_dir_file)
 
     #Use the date.time of the bath file as a basis for
@@ -3725 +3872 @@
     vcur_files = os.listdir(vcur_dir)
     elevation_files.sort()
+
     # the first elevation file should be the
     # file with the same base name as the bath data
@@ -3731 +3879 @@
     number_of_latitudes = bath_grid.shape[0]
     number_of_longitudes = bath_grid.shape[1]
-    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
-
-    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
+
+    longitudes = [bath_metadata['xllcorner'] + x*bath_metadata['cellsize'] \
                   for x in range(number_of_longitudes)]
-    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
+    latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \
                  for y in range(number_of_latitudes)]
 
-     # reverse order of lat, so the fist lat represents the first grid row
+     # reverse order of lat, so the first lat represents the first grid row
     latitudes.reverse()
 
     kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:],
-                                                 minlat=minlat, maxlat=maxlat,
-                                                 minlon=minlon, maxlon=maxlon)
-
+                                                  minlat=minlat, maxlat=maxlat,
+                                                  minlon=minlon, maxlon=maxlon)
 
     bath_grid = bath_grid[kmin:kmax,lmin:lmax]
@@ -3752 +3899 @@
     number_of_longitudes = len(longitudes)
     number_of_times = len(os.listdir(elevation_dir))
-    number_of_points = number_of_latitudes*number_of_longitudes
-    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
+    number_of_points = number_of_latitudes * number_of_longitudes
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
 
     #Work out the times
     if len(elevation_files) > 1:
         # Assume: The time period is less than 24hrs.
-        time_period = (int(elevation_files[1][-3:]) - \
-                      int(elevation_files[0][-3:]))*60*60
+        time_period = (int(elevation_files[1][-3:]) \
+                       - int(elevation_files[0][-3:])) * 60*60
         times = [x*time_period for x in range(len(elevation_files))]
     else:
         times = [0.0]
-
 
     if verbose:
@@ -3769 +3915 @@
         print 'Statistics:'
         print '  Extent (lat/lon):'
-        print '    lat in [%f, %f], len(lat) == %d'\
-              %(min(latitudes), max(latitudes),
-                len(latitudes))
-        print '    lon in [%f, %f], len(lon) == %d'\
-              %(min(longitudes), max(longitudes),
-                len(longitudes))
-        print '    t in [%f, %f], len(t) == %d'\
-              %(min(times), max(times), len(times))
+        print '    lat in [%f, %f], len(lat) == %d' \
+              % (min(latitudes), max(latitudes), len(latitudes))
+        print '    lon in [%f, %f], len(lon) == %d' \
+              % (min(longitudes), max(longitudes), len(longitudes))
+        print '    t in [%f, %f], len(t) == %d' \
+              % (min(times), max(times), len(times))
 
     ######### WRITE THE SWW FILE #############
+
     # NetCDF file definition
     outfile = NetCDFFile(sww_file, 'w')
@@ -3786 +3931 @@
     outfile.description = 'Converted from XXX'
 
-
     #For sww compatibility
     outfile.smoothing = 'Yes'
@@ -3794 +3938 @@
     outfile.starttime = starttime = times[0]
 
-
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
-
     outfile.createDimension('number_of_vertices', 3)
     outfile.createDimension('number_of_points', number_of_points)
@@ -3814 +3956 @@
                                             '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'))
+    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 anuga.coordinate_transforms.redfearn import redfearn
+
     x = zeros(number_of_points, Float)  #Easting
     y = zeros(number_of_points, Float)  #Northing
 
     if verbose: print 'Making triangular grid'
+
     #Get zone of 1st point.
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+    refzone, _, _ = redfearn(latitudes[0], longitudes[0])
 
     vertices = {}
@@ -3842 +3982 @@
     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
+            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)
@@ -3853 +3992 @@
             y[i] = northing
             i += 1
-
 
     #Construct 2 triangles per 'rectangular' element
@@ -3871 +4009 @@
     volumes = array(volumes)
 
-    geo_ref = Geo_reference(refzone,min(x),min(y))
+    geo_ref = Geo_reference(refzone, min(x), min(y))
     geo_ref.write_NetCDF(outfile)
 
     # This will put the geo ref in the middle
-    #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.)
-
+    #geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.)
 
     if verbose:
@@ -3882 +4019 @@
         print 'More Statistics:'
         print '  Extent (/lon):'
-        print '    x in [%f, %f], len(lat) == %d'\
-              %(min(x), max(x),
-                len(x))
-        print '    y in [%f, %f], len(lon) == %d'\
-              %(min(y), max(y),
-                len(y))
-        print 'geo_ref: ',geo_ref
+        print '    x in [%f, %f], len(lat) == %d' \
+              % (min(x), max(x), len(x))
+        print '    y in [%f, %f], len(lon) == %d' \
+              % (min(y), max(y), len(y))
+        print 'geo_ref: ', geo_ref
 
     z = resize(bath_grid,outfile.variables['z'][:].shape)
     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
-    outfile.variables['z'][:] = z # FIXME (Ole): Remove once viewer has been recompiled and changed to use elevation instead of z
-    outfile.variables['elevation'][:] = z  
+# FIXME (Ole): Remove once viewer has been recompiled and changed
+#              to use elevation instead of z
+    outfile.variables['z'][:] = z
+    outfile.variables['elevation'][:] = z
     outfile.variables['volumes'][:] = volumes.astype(Int32) # On Opteron 64
 
@@ -3904 +4041 @@
 
     if verbose: print 'Converting quantities'
+
     n = number_of_times
     for j in range(number_of_times):
         # load in files
         elevation_meta, elevation_grid = \
-           _read_asc(elevation_dir + os.sep + elevation_files[j])
-
-        _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
-        _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
+            _read_asc(elevation_dir + os.sep + elevation_files[j])
+
+        _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j])
+        _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j])
 
         #cut matrix to desired size
@@ -3917 +4055 @@
         u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
         v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
-        
+
         # handle missing values
         missing = (elevation_grid == elevation_meta['NODATA_value'])
         if sometrue (missing):
             if fail_on_NaN:
-                msg = 'File %s contains missing values'\
-                      %(elevation_files[j])
+                msg = 'File %s contains missing values' \
+                      % (elevation_files[j])
                 raise DataMissingValuesError, msg
             else:
-                elevation_grid = elevation_grid*(missing==0) + \
-                                 missing*elevation_NaN_filler
-
-
-        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
+                elevation_grid = elevation_grid*(missing==0) \
+                                 + missing*elevation_NaN_filler
+
+        if verbose and j % ((n+10)/10) == 0: print '  Doing %d of %d' % (j, n)
+
         i = 0
         for k in range(number_of_latitudes):      #Y direction
@@ -3940 +4078 @@
                 ymomentum[j,i] = v_momentum_grid[k,l]*h
                 i += 1
+
     outfile.close()
 
+
+##
+# @brief Return max&min indexes (for slicing) of area specified.
+# @param latitudes_ref ??
+# @param longitudes_ref ??
+# @param minlat Minimum latitude of specified area.
+# @param maxlat Maximum latitude of specified area.
+# @param minlon Minimum longitude of specified area.
+# @param maxlon Maximum longitude of specified area.
+# @return Tuple (lat_min_index, lat_max_index, lon_min_index, lon_max_index)
 def _get_min_max_indexes(latitudes_ref,longitudes_ref,
-                        minlat=None, maxlat=None,
-                        minlon=None, maxlon=None):
-    """
-    return max, min indexes (for slicing) of the lat and long arrays to cover the area
-    specified with min/max lat/long
+                         minlat=None, maxlat=None,
+                         minlon=None, maxlon=None):
+    """
+    Return max, min indexes (for slicing) of the lat and long arrays to cover
+    the area specified with min/max lat/long.
 
     Think of the latitudes and longitudes describing a 2d surface.
@@ -3958 +4107 @@
     assert latitudes are sorted, ascending or decending
     """
+
     latitudes = latitudes_ref[:]
     longitudes = longitudes_ref[:]
@@ -3969 +4119 @@
     #print len(latitudes),len(longitudes)
     #print latitudes[len(latitudes)-1],longitudes[len(longitudes)-1]
-    
+
     lat_ascending = True
     if not allclose(sort(latitudes), latitudes):
         lat_ascending = False
-        # reverse order of lat, so it's in ascending order          
+        # reverse order of lat, so it's in ascending order
         latitudes = latitudes[::-1]
         assert allclose(sort(latitudes), latitudes)
-    #print "latitudes  in funct", latitudes
-    
+
     largest_lat_index = len(latitudes)-1
+
     #Cut out a smaller extent.
     if minlat == None:
@@ -3987 +4137 @@
             lat_min_index = 0
 
-
     if maxlat == None:
         lat_max_index = largest_lat_index #len(latitudes)
@@ -4009 +4158 @@
     # Reversing the indexes, if the lat array is decending
     if lat_ascending is False:
-        lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
+        lat_min_index, lat_max_index = largest_lat_index - lat_max_index, \
                                        largest_lat_index - lat_min_index
     lat_max_index = lat_max_index + 1 # taking into account how slicing works
@@ -4017 +4166 @@
 
 
+##
+# @brief Read an ASC file.
+# @parem filename Path to the file to read.
+# @param verbose True if this function is to be verbose.
 def _read_asc(filename, verbose=False):
     """Read esri file from the following ASCII format (.asc)
@@ -4029 +4182 @@
     NODATA_value  -9999
     138.3698 137.4194 136.5062 135.5558 ..........
-
     """
 
     datafile = open(filename)
 
-    if verbose: print 'Reading DEM from %s' %(filename)
+    if verbose: print 'Reading DEM from %s' % filename
+
     lines = datafile.readlines()
     datafile.close()
@@ -4065 +4218 @@
 
 
-
     ####  URS 2 SWW  ###
 
+# Definitions of various NetCDF dimension names, etc.
 lon_name = 'LON'
 lat_name = 'LAT'
 time_name = 'TIME'
-precision = Float # So if we want to change the precision its done here        
+precision = Float # So if we want to change the precision its done here
+
+##
+# @brief Clas for a NetCDF data file writer.
 class Write_nc:
-    """
-    Write an nc file.
-    
+    """Write an nc file.
+
     Note, this should be checked to meet cdc netcdf conventions for gridded
     data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
-    
-    """
+    """
+
+    ##
+    # @brief Instantiate a Write_nc instance.
+    # @param quantity_name 
+    # @param file_name 
+    # @param time_step_count The number of time steps.
+    # @param time_step The time_step size.
+    # @param lon 
+    # @param lat 
     def __init__(self,
                  quantity_name,
@@ -4087 +4250 @@
                  lon,
                  lat):
-        """
+        """Instantiate a Write_nc instance (NetCDF file writer).
+
         time_step_count is the number of time steps.
         time_step is the time step size
-        
-        pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
+
+        pre-condition: quantity_name must be 'HA', 'UA'or 'VA'.
         """
+
         self.quantity_name = quantity_name
         quantity_units = {'HA':'CENTIMETERS',
-                              'UA':'CENTIMETERS/SECOND',
-                              'VA':'CENTIMETERS/SECOND'}       
-        
-        multiplier_dic = {'HA':100.0, # To convert from m to cm
-                              'UA':100.0,  #  m/s to cm/sec
-                              'VA':-100.0}  # MUX files have positve x in the
-        # Southern direction.  This corrects for it, when writing nc files.
-        
+                          'UA':'CENTIMETERS/SECOND',
+                          'VA':'CENTIMETERS/SECOND'}
+
+        multiplier_dic = {'HA':100.0,   # To convert from m to cm
+                          'UA':100.0,   #             and m/s to cm/sec
+                          'VA':-100.0}  # MUX files have positive x in the
+                                        # Southern direction.  This corrects
+                                        # for it, when writing nc files.
+
         self.quantity_multiplier =  multiplier_dic[self.quantity_name]
-        
+
         #self.file_name = file_name
         self.time_step_count = time_step_count
@@ -4111 +4277 @@
         # NetCDF file definition
         self.outfile = NetCDFFile(file_name, 'w')
-        outfile = self.outfile       
+        outfile = self.outfile
 
         #Create new file
         nc_lon_lat_header(outfile, lon, lat)
-    
+
         # TIME
         outfile.createDimension(time_name, None)
@@ -4123 +4289 @@
         outfile.createVariable(self.quantity_name, precision,
                                (time_name, lat_name, lon_name))
-        outfile.variables[self.quantity_name].missing_value=-1.e+034
-        outfile.variables[self.quantity_name].units= \
+        outfile.variables[self.quantity_name].missing_value = -1.e+034
+        outfile.variables[self.quantity_name].units = \
                                  quantity_units[self.quantity_name]
         outfile.variables[lon_name][:]= ensure_numeric(lon)
@@ -4131 +4297 @@
         #Assume no one will be wanting to read this, while we are writing
         #outfile.close()
-        
+
+    ##
+    # @brief Write a time-step of quantity data.
+    # @param quantity_slice The data to be stored for this time-step.
     def store_timestep(self, quantity_slice):
-        """
-        Write a time slice of quantity info 
+        """Write a time slice of quantity info
+
         quantity_slice is the data to be stored at this time step
         """
-        
-        outfile = self.outfile
-        
+
         # Get the variables
-        time = outfile.variables[time_name]
-        quantity = outfile.variables[self.quantity_name]
-            
+        time = self.outfile.variables[time_name]
+        quantity = self.outfile.variables[self.quantity_name]
+
+        # get index oflice to write
         i = len(time)
 
         #Store time
-        time[i] = i*self.time_step #self.domain.time
-        quantity[i,:] = quantity_slice* self.quantity_multiplier
-        
+        time[i] = i * self.time_step    #self.domain.time
+        quantity[i,:] = quantity_slice * self.quantity_multiplier
+
+    ##
+    # @brief Close file underlying the class instance.
     def close(self):
         self.outfile.close()
 
+
+##
+# @brief Convert URS file to SWW file.
+# @param basename_in Stem of the input filename.
+# @param basename_out Stem of the output filename.
+# @param verbose True if this function is to be verbose.
+# @param remove_nc_files 
+# @param minlat Sets extent of area to be used.  If not supplied, full extent.
+# @param maxlat Sets extent of area to be used.  If not supplied, full extent.
+# @param minlon Sets extent of area to be used.  If not supplied, full extent.
+# @param maxlon Sets extent of area to be used.  If not supplied, full extent.
+# @param mint 
+# @param maxt 
+# @param mean_stage 
+# @param origin A 3-tuple with geo referenced UTM coordinates
+# @param zscale 
+# @param fail_on_NaN 
+# @param NaN_filler 
+# @param elevation 
+# @note Also convert latitude and longitude to UTM. All coordinates are
+#       assumed to be given in the GDA94 datum.
 def urs2sww(basename_in='o', basename_out=None, verbose=False,
             remove_nc_files=True,
             minlat=None, maxlat=None,
-            minlon= None, maxlon=None,
+            minlon=None, maxlon=None,
             mint=None, maxt=None,
             mean_stage=0,
-            origin = None,
+            origin=None,
             zscale=1,
             fail_on_NaN=True,
             NaN_filler=0,
             elevation=None):
-    """
+    """Convert a URS file to an SWW file.
     Convert URS C binary format for wave propagation to
     sww format native to abstract_2d_finite_volumes.
@@ -4179 +4370 @@
     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 assumed to be in decimal degrees. 
+    Lat and lon are assumed to be in decimal degrees.
     NOTE: minlon is the most east boundary.
-    
+
     origin is a 3-tuple with geo referenced
     UTM coordinates (zone, easting, northing)
@@ -4187 +4378 @@
     since all of anuga should be able to handle an arbitary origin.
 
-
     URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
     which means that latitude is the fastest
@@ -4194 +4384 @@
     In URS C binary the latitudes and longitudes are in assending order.
     """
+
     if basename_out == None:
         basename_out = basename_in
+
     files_out = urs2nc(basename_in, basename_out)
+
     ferret2sww(basename_out,
                minlat=minlat,
@@ -4211 +4404 @@
                inverted_bathymetry=True,
                verbose=verbose)
-    #print "files_out",files_out
+    
     if remove_nc_files:
         for file_out in files_out:
             os.remove(file_out)
-    
-def urs2nc(basename_in = 'o', basename_out = 'urs'):
-    """
-    Convert the 3 urs files to 4 nc files.
+
+
+##
+# @brief Convert 3 URS files back to 4 NC files.
+# @param basename_in Stem of the input filenames.
+# @param basename_outStem of the output filenames.
+# @note The name of the urs file names must be:
+#          [basename_in]-z-mux
+#          [basename_in]-e-mux
+#          [basename_in]-n-mux
+def urs2nc(basename_in='o', basename_out='urs'):
+    """Convert the 3 urs files to 4 nc files.
 
     The name of the urs file names must be;
@@ -4224 +4425 @@
     [basename_in]-e-mux
     [basename_in]-n-mux
-    
-    """
-    
+    """
+
     files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
                 basename_in + EAST_VELOCITY_LABEL,
                 basename_in + NORTH_VELOCITY_LABEL]
-    files_out = [basename_out+'_ha.nc',
-                 basename_out+'_ua.nc',
-                 basename_out+'_va.nc']
-    quantities = ['HA','UA','VA']
+    files_out = [basename_out + '_ha.nc',
+                 basename_out + '_ua.nc',
+                 basename_out + '_va.nc']
+    quantities = ['HA', 'UA', 'VA']
 
     #if os.access(files_in[0]+'.mux', os.F_OK) == 0 :
     for i, file_name in enumerate(files_in):
         if os.access(file_name, os.F_OK) == 0:
-            if os.access(file_name+'.mux', os.F_OK) == 0 :
-                msg = 'File %s does not exist or is not accessible' %file_name
+            if os.access(file_name + '.mux', os.F_OK) == 0 :
+                msg = 'File %s does not exist or is not accessible' % file_name
                 raise IOError, msg
             else:
                files_in[i] += '.mux'
                print "file_name", file_name
+
     hashed_elevation = None
     for file_in, file_out, quantity in map(None, files_in,
@@ -4249 +4450 @@
                                            quantities):
         lonlatdep, lon, lat, depth = _binary_c2nc(file_in,
-                                         file_out,
-                                         quantity)
-        #print "lonlatdep", lonlatdep 
+                                                  file_out,
+                                                  quantity)
         if hashed_elevation == None:
-            elevation_file = basename_out+'_e.nc'
+            elevation_file = basename_out + '_e.nc'
             write_elevation_nc(elevation_file,
-                                lon,
-                                lat,
-                                depth)
+                               lon,
+                               lat,
+                               depth)
             hashed_elevation = myhash(lonlatdep)
         else:
             msg = "The elevation information in the mux files is inconsistent"
             assert hashed_elevation == myhash(lonlatdep), msg
+
     files_out.append(elevation_file)
+
     return files_out
-    
+
+
+##
+# @brief Convert a quantity URS file to a NetCDF file.
+# @param file_in Path to input URS file.
+# @param file_out Path to the output file.
+# @param quantity Name of the quantity to be written to the output file.
+# @return A tuple (lonlatdep, lon, lat, depth).
 def _binary_c2nc(file_in, file_out, quantity):
-    """
-    Reads in a quantity urs file and writes a quantity nc file.
-    additionally, returns the depth and lat, long info,
+    """Reads in a quantity urs file and writes a quantity nc file.
+    Additionally, returns the depth and lat, long info,
     so it can be written to a file.
     """
-    columns = 3 # long, lat , depth
+
+    columns = 3                           # long, lat , depth
     mux_file = open(file_in, 'rb')
-    
+
     # Number of points/stations
-    (points_num,)= unpack('i',mux_file.read(4))
+    (points_num,) = unpack('i', mux_file.read(4))
 
     # nt, int - Number of time steps
-    (time_step_count,)= unpack('i',mux_file.read(4))
+    (time_step_count,) = unpack('i', mux_file.read(4))
 
     #dt, float - time step, seconds
     (time_step,) = unpack('f', mux_file.read(4))
-    
+
     msg = "Bad data in the mux file."
     if points_num < 0:
@@ -4293 +4502 @@
         mux_file.close()
         raise ANUGAError, msg
-    
+
     lonlatdep = p_array.array('f')
     lonlatdep.read(mux_file, columns * points_num)
-    lonlatdep = array(lonlatdep, typecode=Float)    
+    lonlatdep = array(lonlatdep, typecode=Float)
     lonlatdep = reshape(lonlatdep, (points_num, columns))
-    
+
     lon, lat, depth = lon_lat2grid(lonlatdep)
     lon_sorted = list(lon)
@@ -4306 +4515 @@
         msg = "Longitudes in mux file are not in ascending order"
         raise IOError, msg
+
     lat_sorted = list(lat)
     lat_sorted.sort()
 
-    if not lat == lat_sorted:
-        msg = "Latitudes in mux file are not in ascending order"
-    
+# UNUSED?
+##    if not lat == lat_sorted:
+##        msg = "Latitudes in mux file are not in ascending order"
+
     nc_file = Write_nc(quantity,
                        file_out,
@@ -4320 +4531 @@
 
     for i in range(time_step_count):
-        #Read in a time slice  from mux file  
+        #Read in a time slice from mux file
         hz_p_array = p_array.array('f')
         hz_p_array.read(mux_file, points_num)
         hz_p = array(hz_p_array, typecode=Float)
         hz_p = reshape(hz_p, (len(lon), len(lat)))
-        hz_p = transpose(hz_p) #mux has lat varying fastest, nc has long v.f. 
+        hz_p = transpose(hz_p)  # mux has lat varying fastest, nc has long v.f.
 
         #write time slice to nc file
         nc_file.store_timestep(hz_p)
+
     mux_file.close()
     nc_file.close()
 
     return lonlatdep, lon, lat, depth
-    
-
+
+
+##
+# @brief Write an NC elevation file.
+# @param file_out Path to the output file.
+# @param lon ??
+# @param lat ??
+# @param depth_vector The elevation data to write.
 def write_elevation_nc(file_out, lon, lat, depth_vector):
-    """
-    Write an nc elevation file.
-    """
-    
+    """Write an nc elevation file."""
+
     # NetCDF file definition
     outfile = NetCDFFile(file_out, 'w')
@@ -4345 +4561 @@
     #Create new file
     nc_lon_lat_header(outfile, lon, lat)
-    
+
     # ELEVATION
     zname = 'ELEVATION'
     outfile.createVariable(zname, precision, (lat_name, lon_name))
-    outfile.variables[zname].units='CENTIMETERS'
-    outfile.variables[zname].missing_value=-1.e+034
-
-    outfile.variables[lon_name][:]= ensure_numeric(lon)
-    outfile.variables[lat_name][:]= ensure_numeric(lat)
-
-    depth = reshape(depth_vector, ( len(lat), len(lon)))
-    outfile.variables[zname][:]= depth
-    
+    outfile.variables[zname].units = 'CENTIMETERS'
+    outfile.variables[zname].missing_value = -1.e+034
+
+    outfile.variables[lon_name][:] = ensure_numeric(lon)
+    outfile.variables[lat_name][:] = ensure_numeric(lat)
+
+    depth = reshape(depth_vector, (len(lat), len(lon)))
+    outfile.variables[zname][:] = depth
+
     outfile.close()
-    
+
+
+##
+# @brief Write lat/lon headers to a NetCDF file.
+# @param outfile Handle to open file to write to.
+# @param lon An iterable of the longitudes.
+# @param lat An iterable of the latitudes.
+# @note Defines lat/long dimensions and variables. Sets various attributes:
+#          .point_spacing  and  .units
+#       and writes lat/lon data.
+
 def nc_lon_lat_header(outfile, lon, lat):
-    """
+    """Write lat/lon headers to a NetCDF file.
+
     outfile is the netcdf file handle.
     lon - a list/array of the longitudes
     lat - a list/array of the latitudes
     """
-    
+
     outfile.institution = 'Geoscience Australia'
     outfile.description = 'Converted from URS binary C'
-    
+
     # Longitude
     outfile.createDimension(lon_name, len(lon))
     outfile.createVariable(lon_name, precision, (lon_name,))
-    outfile.variables[lon_name].point_spacing='uneven'
-    outfile.variables[lon_name].units='degrees_east'
+    outfile.variables[lon_name].point_spacing = 'uneven'
+    outfile.variables[lon_name].units = 'degrees_east'
     outfile.variables[lon_name].assignValue(lon)
-
 
     # Latitude
     outfile.createDimension(lat_name, len(lat))
     outfile.createVariable(lat_name, precision, (lat_name,))
-    outfile.variables[lat_name].point_spacing='uneven'
-    outfile.variables[lat_name].units='degrees_north'
+    outfile.variables[lat_name].point_spacing = 'uneven'
+    outfile.variables[lat_name].units = 'degrees_north'
     outfile.variables[lat_name].assignValue(lat)
 
 
-    
+##
+# @brief 
+# @param long_lat_dep 
+# @return A tuple (long, lat, quantity).
+# @note The latitude is the fastest varying dimension - in mux files.
 def lon_lat2grid(long_lat_dep):
     """
@@ -4397 +4627 @@
     The latitude is the fastest varying dimension - in mux files
     """
+
     LONG = 0
     LAT = 1
@@ -4402 +4633 @@
 
     long_lat_dep = ensure_numeric(long_lat_dep, Float)
-    
+
     num_points = long_lat_dep.shape[0]
     this_rows_long = long_lat_dep[0,LONG]
@@ -4410 +4641 @@
     while long_lat_dep[i,LONG] == this_rows_long and i < num_points:
         i += 1
+
     # determine the lats and longsfrom the grid
-    lat = long_lat_dep[:i, LAT]        
+    lat = long_lat_dep[:i, LAT]
     long = long_lat_dep[::i, LONG]
-    
+
     lenlong = len(long)
     lenlat = len(lat)
-    #print 'len lat', lat, len(lat)
-    #print 'len long', long, len(long) 
-          
-    msg = 'Input data is not gridded'      
+
+    msg = 'Input data is not gridded'
     assert num_points % lenlat == 0, msg
     assert num_points % lenlong == 0, msg
-          
-    # Test that data is gridded        
+
+    # Test that data is gridded
     for i in range(lenlong):
         msg = 'Data is not gridded.  It must be for this operation'
-        first = i*lenlat
+        first = i * lenlat
         last = first + lenlat
-                
+
         assert allclose(long_lat_dep[first:last,LAT], lat), msg
         assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
-    
-    
-#    print 'range long', min(long), max(long)
-#    print 'range lat', min(lat), max(lat)
-#    print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0])
-#    print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1])
-    
-   
-    
+
     msg = 'Out of range latitudes/longitudes'
     for l in lat:assert -90 < l < 90 , msg
@@ -4448 +4670 @@
     # FIXME - make this faster/do this a better way
     # use numeric transpose, after reshaping the quantity vector
-#    quantity = zeros(len(long_lat_dep), Float)
     quantity = zeros(num_points, Float)
-    
-#    print 'num',num_points
+
     for lat_i, _ in enumerate(lat):
         for long_i, _ in enumerate(long):
-            q_index = lat_i*lenlong+long_i
-            lld_index = long_i*lenlat+lat_i
-#            print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index
+            q_index = lat_i*lenlong + long_i
+            lld_index = long_i*lenlat + lat_i
             temp = long_lat_dep[lld_index, QUANTITY]
             quantity[q_index] = temp
-            
+
     return long, lat, quantity
 
-####  END URS 2 SWW  ###
-
-#### URS UNGRIDDED 2 SWW ###
+################################################################################
+# END URS 2 SWW
+################################################################################
+
+################################################################################
+# URS UNGRIDDED 2 SWW
+################################################################################
 
 ### PRODUCING THE POINTS NEEDED FILE ###
@@ -4475 +4698 @@
 #LONG_AMOUNT = 3600
 
+
+##
+# @brief 
+# @param file_name 
+# @param boundary_polygon 
+# @param zone 
+# @param ll_lat 
+# @param ll_long 
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param isSouthernHemisphere 
+# @param export_csv 
+# @param use_cache 
+# @param verbose True if this function is to be verbose.
+# @return 
 def URS_points_needed_to_file(file_name, boundary_polygon, zone,
                               ll_lat, ll_long,
-                              grid_spacing, 
+                              grid_spacing,
                               lat_amount, long_amount,
                               isSouthernHemisphere=True,
@@ -4489 +4728 @@
     1st line is the number of points,
     each line after represents a point,in lats and longs.
-    
+
     Note: The polygon cannot cross zones or hemispheres.
 
     A work-a-round for different zones or hemispheres is to run this twice,
     once for each zone, and then combine the output.
-    
+
     file_name - name of the urs file produced for David.
     boundary_polygon - a list of points that describes a polygon.
@@ -4501 +4740 @@
 
      This is info about the URS model that needs to be inputted.
-     
+
     ll_lat - lower left latitude, in decimal degrees
     ll-long - lower left longitude, in decimal degrees
     grid_spacing - in deciamal degrees
-    lat_amount - number of latitudes 
-    long_amount- number of longs 
-
+    lat_amount - number of latitudes
+    long_amount- number of longs
 
     Don't add the file extension.  It will be added.
     """
+
     geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long,
-                            grid_spacing, 
+                            grid_spacing,
                             lat_amount, long_amount, isSouthernHemisphere,
                             use_cache, verbose)
+
     if not file_name[-4:] == ".urs":
         file_name += ".urs"
+
     geo.export_points_file(file_name, isSouthHemisphere=isSouthernHemisphere)
+
     if export_csv:
         if file_name[-4:] == ".urs":
             file_name = file_name[:-4] + ".csv"
         geo.export_points_file(file_name)
+
     return geo
 
+
+##
+# @brief 
+# @param boundary_polygon
+# @param zone
+# @param ll_lat
+# @param ll_long
+# @param grid_spacing
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
+# @param use_cache
+# @param verbose
 def URS_points_needed(boundary_polygon, zone, ll_lat,
-                      ll_long, grid_spacing, 
+                      ll_long, grid_spacing,
                       lat_amount, long_amount, isSouthHemisphere=True,
                       use_cache=False, verbose=False):
     args = (boundary_polygon,
             zone, ll_lat,
-            ll_long, grid_spacing, 
+            ll_long, grid_spacing,
             lat_amount, long_amount, isSouthHemisphere)
-    kwargs = {}  
+    kwargs = {}
+
     if use_cache is True:
         try:
             from anuga.caching import cache
         except:
-            msg = 'Caching was requested, but caching module'+\
+            msg = 'Caching was requested, but caching module' \
                   'could not be imported'
             raise msg
 
-
         geo = cache(_URS_points_needed,
-                  args, kwargs,
-                  verbose=verbose,
-                  compression=False)
+                    args, kwargs,
+                    verbose=verbose,
+                    compression=False)
     else:
         geo = apply(_URS_points_needed, args, kwargs)
@@ -4551 +4807 @@
     return geo
 
+
+##
+# @brief 
+# @param boundary_polygon An iterable of points that describe a polygon.
+# @param zone
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing Grid spacing in decimal degrees.
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
 def _URS_points_needed(boundary_polygon,
-                      zone, ll_lat,
-                      ll_long, grid_spacing, 
-                      lat_amount, long_amount,
+                       zone, ll_lat,
+                       ll_long, grid_spacing,
+                       lat_amount, long_amount,
                        isSouthHemisphere):
     """
@@ -4563 +4830 @@
     ll_lat - lower left latitude, in decimal degrees
     ll-long - lower left longitude, in decimal degrees
-    grid_spacing - in deciamal degrees
-
-    """
-    
+    grid_spacing - in decimal degrees
+
+    """
+
     from sets import ImmutableSet
-    
+
     msg = "grid_spacing can not be zero"
-    assert not grid_spacing == 0, msg 
+    assert not grid_spacing == 0, msg
+
     a = boundary_polygon
+
     # List of segments.  Each segment is two points.
     segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
+
     # convert the segs to Lat's and longs.
-    
     # Don't assume the zone of the segments is the same as the lower left
     # corner of the lat long data!!  They can easily be in different zones
-    
     lat_long_set = ImmutableSet()
     for seg in segs:
-        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, 
-                      lat_amount, long_amount, zone, isSouthHemisphere)
+        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
+                                        lat_amount, long_amount, zone,
+                                        isSouthHemisphere)
         lat_long_set |= ImmutableSet(points_lat_long)
+
     if lat_long_set == ImmutableSet([]):
-        msg = """URS region specified and polygon does not overlap."""
+        msg = "URS region specified and polygon does not overlap."
         raise ValueError, msg
 
@@ -4591 +4861 @@
     # these points.  There should be.
     geo = Geospatial_data(data_points=list(lat_long_set),
-                              points_are_lats_longs=True)
+                          points_are_lats_longs=True)
+
     return geo
-    
-def points_needed(seg, ll_lat, ll_long, grid_spacing, 
+
+
+##
+# @brief Get the points that are needed to interpolate any point a a segment.
+# @param seg Two points in the UTM.
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param zone 
+# @param isSouthHemisphere 
+# @return A list of points.
+def points_needed(seg, ll_lat, ll_long, grid_spacing,
                   lat_amount, long_amount, zone,
                   isSouthHemisphere):
@@ -4602 +4885 @@
     interpolate any point on the segment.
     """
+
     from math import sqrt
-    #print "zone",zone 
+
     geo_reference = Geo_reference(zone=zone)
-    #print "seg", seg
-    geo = Geospatial_data(seg,geo_reference=geo_reference)
+    geo = Geospatial_data(seg, geo_reference=geo_reference)
     seg_lat_long = geo.get_data_points(as_lat_long=True,
                                        isSouthHemisphere=isSouthHemisphere)
-    #print "seg_lat_long", seg_lat_long
+
     # 1.415 = 2^0.5, rounded up....
     sqrt_2_rounded_up = 1.415
     buffer = sqrt_2_rounded_up * grid_spacing
-    
+
     max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
     max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
@@ -4619 +4902 @@
     min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
 
-    #print "min_long", min_long
-    #print "ll_long", ll_long
-    #print "grid_spacing", grid_spacing
-    first_row = (min_long - ll_long)/grid_spacing
+    first_row = (min_long - ll_long) / grid_spacing
+
     # To round up
     first_row_long = int(round(first_row + 0.5))
-    #print "first_row", first_row_long
-
-    last_row = (max_long - ll_long)/grid_spacing # round down
+
+    last_row = (max_long - ll_long) / grid_spacing # round down
     last_row_long = int(round(last_row))
-    #print "last_row",last_row _long
-    
-    first_row = (min_lat - ll_lat)/grid_spacing
+
+    first_row = (min_lat - ll_lat) / grid_spacing
     # To round up
     first_row_lat = int(round(first_row + 0.5))
-    #print "first_row", first_row_lat
-
-    last_row = (max_lat - ll_lat)/grid_spacing # round down
+
+    last_row = (max_lat - ll_lat) / grid_spacing # round down
     last_row_lat = int(round(last_row))
-    #print "last_row",last_row_lat
-
-    # to work out the max distance -
-    # 111120 - horizontal distance between 1 deg latitude. 
-    #max_distance = sqrt_2_rounded_up * 111120 * grid_spacing
+
     max_distance = 157147.4112 * grid_spacing
-    #print "max_distance", max_distance #2619.12 m for 1 minute
     points_lat_long = []
+
     # Create a list of the lat long points to include.
     for index_lat in range(first_row_lat, last_row_lat + 1):
         for index_long in range(first_row_long, last_row_long + 1):
-            #print "index_lat", index_lat
-            #print "index_long", index_long
             lat = ll_lat + index_lat*grid_spacing
             long = ll_long + index_long*grid_spacing
@@ -4657 +4929 @@
             if keep_point(lat, long, seg, max_distance):
                 points_lat_long.append((lat, long)) #must be hashable
-    #print "points_lat_long", points_lat_long
 
     # Now that we have these points, lets throw ones out that are too far away
     return points_lat_long
 
+
+##
+# @brief 
+# @param lat
+# @param long
+# @param seg Two points in UTM.
+# @param max_distance
 def keep_point(lat, long, seg, max_distance):
     """
     seg is two points, UTM
     """
+
     from math import sqrt
+
     _ , x0, y0 = redfearn(lat, long)
     x1 = seg[0][0]
@@ -4678 +4958 @@
             (x2_1)*(x2_1)+(y2_1)*(y2_1))
     except ZeroDivisionError:
-        #print "seg", seg
-        #print "x0", x0
-        #print "y0", y0
-        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 and \
-           abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
+        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 \
+           and abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
             return True
         else:
             return False
-    
-    if d <= max_distance:
-        return True
-    else:
-        return False
-    
-    #### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE ####
-   
+
+    return d <= max_distance
+
+################################################################################
+# CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
+################################################################################
+
 WAVEHEIGHT_MUX_LABEL = '-z-mux'
 EAST_VELOCITY_LABEL =  '-e-mux'
-NORTH_VELOCITY_LABEL =  '-n-mux' 
+NORTH_VELOCITY_LABEL =  '-n-mux'
+
+##
+# @brief Convert URS file(s) (wave prop) to an SWW file.
+# @param basename_in Stem of the input filenames.
+# @param basename_out Path to the output SWW file.
+# @param verbose True if this function is to be verbose.
+# @param mint
+# @param maxt
+# @param mean_stage
+# @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing).
+# @param hole_points_UTM
+# @param zscale
+# @note Also convert latitude and longitude to UTM. All coordinates are
+#       assumed to be given in the GDA94 datum.
+# @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux'
+#       added for relative height, x-velocity and y-velocity respectively.
 def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
                       mint=None, maxt=None,
@@ -4703 +4995 @@
                       hole_points_UTM=None,
                       zscale=1):
-    """   
+    """
     Convert URS C binary format for wave propagation to
     sww format native to abstract_2d_finite_volumes.
-
 
     Specify only basename_in and read files of the form
@@ -4719 +5010 @@
     min's and max's: If omitted - full extend is used.
     To include a value min ans max may equal it.
-    Lat and lon are assumed to be in decimal degrees. 
-    
+    Lat and lon are assumed to be in decimal degrees.
+
     origin is a 3-tuple with geo referenced
     UTM coordinates (zone, easting, northing)
@@ -4727 +5018 @@
     The mux point info is NOT relative to this origin.
 
-
-    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
+    URS C binary format has data organised as TIME, LONGITUDE, LATITUDE
     which means that latitude is the fastest
     varying dimension (row major order, so to speak)
@@ -4738 +5028 @@
     function used, and the different grid structure between urs2sww
     and this function.
-    
+
     Interpolating data that has an underlying gridded source can
     easily end up with different values, depending on the underlying
@@ -4750 +5040 @@
     The grid can be
      -
-    |\|    A
+    |\|   A
      -
      or;
-      - 
-     |/|   B 
       -
-      If a point is just below the center of the midpoint, it will have a
-      +ve value in grid A and a -ve value in grid B.
-    """ 
+     |/|  B
+      -
+
+    If a point is just below the center of the midpoint, it will have a
+    +ve value in grid A and a -ve value in grid B.
+    """
+
     from anuga.mesh_engine.mesh_engine import NoTrianglesError
     from anuga.pmesh.mesh import Mesh
@@ -4767 +5059 @@
     quantities = ['HA','UA','VA']
 
-    # instanciate urs_points of the three mux files.
+    # instantiate urs_points of the three mux files.
     mux = {}
     for quantity, file in map(None, quantities, files_in):
         mux[quantity] = Urs_points(file)
-        
+
     # Could check that the depth is the same. (hashing)
 
     # handle to a mux file to do depth stuff
     a_mux = mux[quantities[0]]
-    
+
     # Convert to utm
     lat = a_mux.lonlatdep[:,1]
     long = a_mux.lonlatdep[:,0]
-    points_utm, zone = convert_from_latlon_to_utm( \
-        latitudes=lat, longitudes=long)
-    #print "points_utm", points_utm
-    #print "zone", zone
-
-    elevation = a_mux.lonlatdep[:,2] * -1 #
-    
-    # grid ( create a mesh from the selected points)
+    points_utm, zone = convert_from_latlon_to_utm(latitudes=lat,
+                                                  longitudes=long)
+
+    elevation = a_mux.lonlatdep[:,2] * -1
+
+    # grid (create a mesh from the selected points)
     # This mesh has a problem.  Triangles are streched over ungridded areas.
-    #  If these areas could be described as holes in pmesh, that would be great
+    # If these areas could be described as holes in pmesh, that would be great.
 
     # I can't just get the user to selection a point in the middle.
@@ -4798 +5088 @@
     mesh.add_vertices(points_utm)
     mesh.auto_segment(smooth_indents=True, expand_pinch=True)
+
     # To try and avoid alpha shape 'hugging' too much
-    mesh.auto_segment( mesh.shape.get_alpha()*1.1 )
+    mesh.auto_segment(mesh.shape.get_alpha() * 1.1)
     if hole_points_UTM is not None:
         point = ensure_absolute(hole_points_UTM)
         mesh.add_hole(point[0], point[1])
+
     try:
         mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
     except NoTrianglesError:
-        # This is a bit of a hack, going in and changing the
-        # data structure.
+        # This is a bit of a hack, going in and changing the data structure.
         mesh.holes = []
         mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
+
     mesh_dic = mesh.Mesh2MeshList()
 
     #mesh.export_mesh_file(basename_in + '_168.tsh')
-    #import sys; sys.exit() 
+    #import sys; sys.exit()
     # These are the times of the mux file
     mux_times = []
     for i in range(a_mux.time_step_count):
-        mux_times.append(a_mux.time_step * i)  
-    mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt)
+        mux_times.append(a_mux.time_step * i)
+    (mux_times_start_i, mux_times_fin_i) = mux2sww_time(mux_times, mint, maxt)
     times = mux_times[mux_times_start_i:mux_times_fin_i]
-    
+
     if mux_times_start_i == mux_times_fin_i:
         # Close the mux files
         for quantity, file in map(None, quantities, files_in):
             mux[quantity].close()
-        msg="Due to mint and maxt there's no time info in the boundary SWW."
+        msg = "Due to mint and maxt there's no time info in the boundary SWW."
         raise Exception, msg
-        
+
     # If this raise is removed there is currently no downstream errors
-           
+
     points_utm=ensure_numeric(points_utm)
-    assert ensure_numeric(mesh_dic['generatedpointlist']) == \
-           ensure_numeric(points_utm)
-    
+    assert ensure_numeric(mesh_dic['generatedpointlist']) \
+           == ensure_numeric(points_utm)
+
     volumes = mesh_dic['generatedtrianglelist']
-    
-    # write sww intro and grid stuff.   
+
+    # write sww intro and grid stuff.
     if basename_out is None:
         swwname = basename_in + '.sww'
@@ -4843 +5135 @@
 
     if verbose: print 'Output to ', swwname
+
     outfile = NetCDFFile(swwname, 'w')
+
     # For a different way of doing this, check out tsh2sww
     # work out sww_times and the index range this covers
     sww = Write_sww()
     sww.store_header(outfile, times, len(volumes), len(points_utm),
-                     verbose=verbose,sww_precision=Float)
+                     verbose=verbose, sww_precision=Float)
     outfile.mean_stage = mean_stage
     outfile.zscale = zscale
@@ -4855 +5149 @@
                             elevation, zone,  new_origin=origin,
                             verbose=verbose)
-    
+
     if verbose: print 'Converting quantities'
+
+    # Read in a time slice from each mux file and write it to the SWW file
     j = 0
-    # Read in a time slice from each mux file and write it to the sww file
     for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
         if j >= mux_times_start_i and j < mux_times_fin_i:
@@ -4864 +5159 @@
             h = stage - elevation
             xmomentum = ua*h
-            ymomentum = -1*va*h # -1 since in mux files south is positive.
-            sww.store_quantities(outfile, 
-                                 slice_index=j - mux_times_start_i,
+            ymomentum = -1 * va * h # -1 since in mux files south is positive.
+            sww.store_quantities(outfile,
+                                 slice_index=j-mux_times_start_i,
                                  verbose=verbose,
                                  stage=stage,
@@ -4873 +5168 @@
                                  sww_precision=Float)
         j += 1
+
     if verbose: sww.verbose_quantities(outfile)
+
     outfile.close()
-    #
+
     # Do some conversions while writing the sww file
 
-    ##################################
-    # READ MUX2 FILES line of points #
-    ##################################
+
+################################################################################
+# READ MUX2 FILES line of points
+################################################################################
 
 WAVEHEIGHT_MUX2_LABEL = '-z-mux2'
-EAST_VELOCITY_MUX2_LABEL =  '-e-mux2'
-NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'    
-
+EAST_VELOCITY_MUX2_LABEL = '-e-mux2'
+NORTH_VELOCITY_MUX2_LABEL = '-n-mux2'
+
+##
+# @brief 
+# @param filenames List of mux2 format input filenames.
+# @param weights Weights associated with each source file.
+# @param permutation The gauge numbers for which data is extracted.
+# @param verbose True if this function is to be verbose.
+# @return (times, latitudes, longitudes, elevation, quantity, starttime)
 def read_mux2_py(filenames,
                  weights=None,
                  permutation=None,
                  verbose=False):
-    """Access the mux files specified in the filenames list. Combine the 
-       data found therin as a weighted linear sum as specifed by the weights. 
-       If permutation is None or empty extract timeseries data for all gauges within the 
-       files. 
-         
-       Input: 
-           filenames:   List of filenames specifiying the file containing the 
-                        timeseries data (mux2 format) for each source 
-           weights:     Weighs associated with each source (defaults to 1 for each source)
-           permutation: Specifies the gauge numbers that for which data is to be 
-                        extracted 
-    """ 
+    """Access the mux files specified in the filenames list. Combine the
+       data found therin as a weighted linear sum as specifed by the weights.
+       If permutation is None or empty extract timeseries data for all gauges
+       within the files.
+
+       Input:
+           filenames:   List of filenames specifiying the file containing the
+                        timeseries data (mux2 format) for each source
+           weights:     Weighs associated with each source
+                        (defaults to 1 for each source)
+           permutation: Specifies the gauge numbers that for which data is to be
+                        extracted
+    """
 
     from Numeric import ones,Float,compress,zeros,arange
     from urs_ext import read_mux2
 
-    numSrc=len(filenames)
-    
-    file_params=-1*ones(3,Float) #[nsta,dt,nt]
-    
+    numSrc = len(filenames)
+
+    file_params = -1 * ones(3, Float)                    # [nsta,dt,nt]
+
     # Convert verbose to int C flag
     if verbose:
@@ -4918 +5224 @@
     if weights is None:
         weights = ones(numSrc)
-        
+
     if permutation is None:
-        permutation = ensure_numeric([], Float)    
-
-    #print 'filenames', filenames
-    #print 'weights', weights
-                
-    # Call underlying C implementation urs2sts_ext.c