Changeset 7735

Timestamp:

May 20, 2010, 12:36:35 PM (15 years ago)

Author:

James Hudson

Message:

Split up some of the huge modules in shallow_water, fixed most of the unit test dependencies.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/test_gauge.py (modified) (2 diffs)
• abstract_2d_finite_volumes/test_generic_boundary_conditions.py (modified) (1 diff)
• abstract_2d_finite_volumes/test_util.py (modified) (2 diffs)
• abstract_2d_finite_volumes/util.py (modified) (2 diffs)
• culvert_flows/culvert_class.py (modified) (1 diff)
• culvert_flows/new_culvert_class.py (modified) (1 diff)
• damage_modelling/test_inundation_damage.py (modified) (2 diffs)
• fit_interpolate/test_interpolate.py (modified) (1 diff)
• interface/__init__.py (deleted)
• interface/boundaries.py (deleted)
• interface/model.py (deleted)
• interface/test_all.py (deleted)
• interface/test_model.py (deleted)
• shallow_water/boundaries.py (modified) (3 diffs)
• shallow_water/data_manager.py (modified) (13 diffs)
• shallow_water/shallow_water_domain.py (modified) (1 diff)
• shallow_water/sww_file.py (added)
• shallow_water/test_data_manager.py (modified) (15 diffs)
• shallow_water/test_forcing.py (modified) (4 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (1 diff)
• shallow_water/tsh2sww.py (modified) (1 diff)
• utilities/file_utils.py (added)

anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

@@ -10 +10 @@
 from anuga.abstract_2d_finite_volumes.gauge import *
 from anuga.config import epsilon
-from anuga.shallow_water.data_manager import timefile2netcdf,del_dir
 
 from anuga.utilities.numerical_tools import NAN,mean
@@ -18 +17 @@
 from anuga.pmesh.mesh import Mesh
 from anuga.shallow_water import Domain, Transmissive_boundary
-from anuga.shallow_water.data_manager import SWW_file
+from anuga.shallow_water.sww_file import SWW_file
+from anuga.shallow_water.file_conversion import timefile2netcdf
+from anuga.utilities.file_utils import del_dir
 from csv import reader,writer
 import time

anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

@@ -404 +404 @@
 
         #Convert ASCII file to NetCDF (Which is what we really like!)
-        from anuga.shallow_water.data_manager import timefile2netcdf
+        from anuga.shallow_water.file_conversion import timefile2netcdf
         
         timefile2netcdf(filename, quantity_names = ['stage', 'xmomentum'])

anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

@@ -9 +9 @@
 from anuga.abstract_2d_finite_volumes.util import *
 from anuga.config import epsilon
-from anuga.shallow_water.data_manager import timefile2netcdf,del_dir
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.shallow_water.file_conversion import timefile2netcdf
+from anuga.utilities.file_utils import del_dir
 
 from anuga.utilities.numerical_tools import NAN
@@ -17 +19 @@
 from anuga.pmesh.mesh import Mesh
 from anuga.shallow_water import Domain, Transmissive_boundary
-from anuga.shallow_water.data_manager import SWW_file
+from anuga.shallow_water.sww_file import SWW_file
 from csv import reader,writer
 import time

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -650 +650 @@
 
     from os import sep
-    from anuga.shallow_water.data_manager import \
-                               get_all_files_with_extension, csv2dict
+    from anuga.utilities.file_utils import load_csv_as_dict, \
+                                get_all_files_with_extension
 
     seconds_in_hour = 3600
@@ -955 +955 @@
     """
 
-    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
-                                                 get_maximum_inundation_data,\
-                                                 csv2dict
+    from anuga.shallow_water.data_manager import \
+        get_maximum_inundation_data, csv2dict
                                                  
     file = open(runup_filename, "w")

anuga_core/source/anuga/culvert_flows/culvert_class.py

@@ -1 +1 @@
 import sys
 
-from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
+from anuga.shallow_water.forcing import Inflow, General_forcing
 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
 from anuga.utilities.system_tools import log_to_file

anuga_core/source/anuga/culvert_flows/new_culvert_class.py

@@ -1 +1 @@
 import sys
 
-from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
+from anuga.shallow_water.forcing import Inflow, General_forcing
 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
 from anuga.utilities.system_tools import log_to_file

anuga_core/source/anuga/damage_modelling/test_inundation_damage.py

@@ -15 +15 @@
 from anuga.shallow_water import Domain, Transmissive_boundary
 from anuga.utilities.numerical_tools import mean
-from anuga.shallow_water.data_manager import SWW_file
+from anuga.shallow_water.sww_file import SWW_file
 
 import numpy as num
@@ -293 +293 @@
         fd = open(csv_file,'wb')
         writer = csv.writer(fd)
-        writer.writerow(['x','y',STR_VALUE_LABEL,CONT_VALUE_LABEL,'ROOF_TYPE',WALL_TYPE_LABEL, SHORE_DIST_LABEL])
+        writer.writerow(['x', 'y', STR_VALUE_LABEL, CONT_VALUE_LABEL, \
+        'ROOF_TYPE', WALL_TYPE_LABEL, SHORE_DIST_LABEL])
         writer.writerow([5.5,0.5,'10','130000','Metal','Timber',20])
         writer.writerow([4.5,1.0,'150','76000','Metal','Double Brick',20])

anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -24 +24 @@
 from anuga.shallow_water import Domain, Transmissive_boundary
 from anuga.utilities.numerical_tools import mean, NAN
-from anuga.shallow_water.data_manager import SWW_file
+from anuga.shallow_water.sww_file import SWW_file
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 from anuga.pmesh.mesh import Mesh

anuga_core/source/anuga/shallow_water/boundaries.py

@@ -188 +188 @@
 
     def __init__(self, domain=None, function=None):
-        """ Instantiate a Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
+        """ Instantiate a
+            Transmissive_n_momentum_zero_t_momentum_set_stage_boundary.
             domain is the domain containing the boundary
             function is the function to apply
@@ -209 +210 @@
     def __repr__(self):
         """ Return a representation of this instance. """
-        return 'Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(%s)' %self.domain
+        msg='Transmissive_n_momentum_zero_t_momentum_set_stage_boundary'
+        msg+='(%s)' %self.domain
+        return msg
 
 
@@ -393 +396 @@
             
         # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
-        # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain. 
+        # Where v is velocity, n is manning's coefficient, h is depth
+        # and S is the slope into the domain. 
         # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S) 
         # from which we can isolate depth to get

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -94 +94 @@
 import anuga.utilities.log as log
 
+from anuga.utilities.file_utils import create_filename,\
+                        get_all_swwfiles, load_csv_as_dict
+from sww_file import Read_sww, Write_sww
+
 
 # Default block size for sww2dem()
@@ -118 +122 @@
                     'speed': \
  '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
-
-
-##
-# @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
-    """
-
-    s = s.strip()
-    s = s.replace(' ', '_')
-    s = s.replace('(', '')
-    s = s.replace(')', '')
-    s = s.replace('__', '_')
-
-    return s
-
-
-##
-# @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.
-    If path does not exist it will be created if possible
-
-    USAGE:
-       checkdir(path, verbose):
-
-    ARGUMENTS:
-        path -- Directory
-        verbose -- Flag verbose output (default: None)
-
-    RETURN VALUE:
-        Verified path including trailing separator
-    """
-
-    import os.path
-
-    if sys.platform in ['nt', 'dos', 'win32', 'what else?']:
-        unix = 0
-    else:
-        unix = 1
-
-    # add terminal separator, if it's not already there
-    if path[-1] != os.sep:
-        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)
-
-            # Change access rights if possible
-            if unix:
-                exitcode = os.system('chmod 775 ' + path)
-            else:
-                pass  # FIXME: What about access rights under Windows?
-
-            if verbose: log.critical('MESSAGE: Directory %s created.' % path)
-        except:
-            log.critical('WARNING: Directory %s could not be created.' % path)
-            if unix:
-                path = '/tmp/'
-            else:
-                path = 'C:' + os.sep
-
-            log.critical("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
-    """
-
-    if os.path.isdir(path):
-        for file in os.listdir(path):
-            X = os.path.join(path, file)
-
-            if os.path.isdir(X) and not os.path.islink(X):
-                del_dir(X)
-            else:
-                try:
-                    os.remove(X)
-                except:
-                    log.critical("Could not remove file %s" % X)
-
-        os.rmdir(path)
-
-
-##
-# @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)
-        if verbose: log.critical('Removed %s' % path)
-    except OSError, (errno, strerror):
-        log.critical(ERROR_STR % {'path' : path, 'error': strerror })
-
-
-##
-# @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
-
-    for x in os.listdir(path):
-        fullpath = os.path.join(path, x)
-        if os.path.isfile(fullpath):
-            f = os.remove
-            rmgeneric(fullpath, f)
-        elif os.path.isdir(fullpath):
-            removeall(fullpath)
-            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):
-    FN = check_dir(datadir) + filename
-
-    if size is not None:
-        FN += '_size%d' % size
-
-    if time is not None:
-        FN += '_time%.2f' % time
-
-    FN += '.' + format
-
-    return FN
-
-
-##
-# @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
-    """
-
-    import glob
-
-    dir = check_dir(datadir)
-    pattern = dir + os.sep + filename + '_size=%d*.%s' % (size, format)
-
-    return glob.glob(pattern)
-
-
-##
-# @brief Generic class for storing output to e.g. visualisation or checkpointing
-class Data_format:
-    """Generic interface to data formats
-    """
-
-    ##
-    # @brief Instantiate this instance.
-    # @param domain 
-    # @param extension 
-    # @param mode The mode of the underlying file.
-    def __init__(self, domain, extension, mode=netcdf_mode_w):
-        assert mode[0] in ['r', 'w', 'a'], \
-               "Mode %s must be either:\n" % mode + \
-               "   'w' (write)\n" + \
-               "   'r' (read)\n" + \
-               "   'a' (append)"
-
-        #Create filename
-        self.filename = create_filename(domain.get_datadir(),
-                                        domain.get_name(), extension)
-
-        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_volumes = domain.number_of_full_triangles
-        #self.number_of_volumes = len(domain)
-
-        #FIXME: Should we have a general set_precision function?
-
-
-##
-# @brief Class for storing output to e.g. visualisation
-class SWW_file(Data_format):
-    """Interface to native NetCDF format (.sww) for storing model output
-
-    There are two kinds of data
-
-    1: Constant data: Vertex coordinates and field values. Stored once
-    2: Variable data: Conserved quantities. Stored once per timestep.
-
-    All data is assumed to reside at vertex locations.
-    """
-
-    ##
-    # @brief Instantiate this instance.
-    # @param domain ??
-    # @param mode Mode of the underlying data file.
-    # @param max_size ??
-    # @param recursion ??
-    # @note Prepare the underlying data file if mode starts with 'w'.
-    def __init__(self, domain, 
-                 mode=netcdf_mode_w, max_size=2000000000, recursion=False):
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        self.precision = netcdf_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
-        if hasattr(domain, 'minimum_storable_height'):
-            self.minimum_storable_height = domain.minimum_storable_height
-        else:
-            self.minimum_storable_height = default_minimum_storable_height
-
-        # Call parent constructor
-        Data_format.__init__(self, domain, 'sww', mode)
-
-        # Get static and dynamic quantities from domain
-        static_quantities = []
-        dynamic_quantities = []
-        
-        for q in domain.quantities_to_be_stored:
-            flag = domain.quantities_to_be_stored[q]
-        
-            msg = 'Quantity %s is requested to be stored ' % q
-            msg += 'but it does not exist in domain.quantities'
-            assert q in domain.quantities, msg
-        
-            assert flag in [1,2]
-            if flag == 1: static_quantities.append(q)
-            if flag == 2: dynamic_quantities.append(q)                
-                       
-        
-        # NetCDF file definition
-        fid = NetCDFFile(self.filename, mode)
-        if mode[0] == 'w':
-            description = 'Output from anuga.abstract_2d_finite_volumes ' \
-                          'suitable for plotting'
-                          
-            self.writer = Write_sww(static_quantities, dynamic_quantities)
-            self.writer.store_header(fid,
-                                     domain.starttime,
-                                     self.number_of_volumes,
-                                     self.domain.number_of_full_nodes,
-                                     description=description,
-                                     smoothing=domain.smooth,
-                                     order=domain.default_order,
-                                     sww_precision=self.precision)
-
-            # Extra optional information
-            if hasattr(domain, 'texture'):
-                fid.texture = domain.texture
-
-            if domain.quantities_to_be_monitored is not None:
-                fid.createDimension('singleton', 1)
-                fid.createDimension('two', 2)
-
-                poly = domain.monitor_polygon
-                if poly is not None:
-                    N = len(poly)
-                    fid.createDimension('polygon_length', N)
-                    fid.createVariable('extrema.polygon',
-                                       self.precision,
-                                       ('polygon_length', 'two'))
-                    fid.variables['extrema.polygon'][:] = poly
-
-                interval = domain.monitor_time_interval
-                if interval is not None:
-                    fid.createVariable('extrema.time_interval',
-                                       self.precision,
-                                       ('two',))
-                    fid.variables['extrema.time_interval'][:] = interval
-
-                for q in domain.quantities_to_be_monitored:
-                    fid.createVariable(q + '.extrema', self.precision,
-                                       ('numbers_in_range',))
-                    fid.createVariable(q + '.min_location', self.precision,
-                                       ('numbers_in_range',))
-                    fid.createVariable(q + '.max_location', self.precision,
-                                       ('numbers_in_range',))
-                    fid.createVariable(q + '.min_time', self.precision,
-                                       ('singleton',))
-                    fid.createVariable(q + '.max_time', self.precision,
-                                       ('singleton',))
-
-        fid.close()
-
-    ##
-    # @brief Store connectivity data into the underlying data file.
-    def store_connectivity(self):
-        """Store information about nodes, triangles and static quantities
-
-        Writes x,y coordinates of triangles and their connectivity.
-        
-        Store also any quantity that has been identified as static.
-        """
-
-        # FIXME: Change name to reflect the fact thta this function 
-        # stores both connectivity (triangulation) and static quantities
-        
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        domain = self.domain
-
-        # append to the NetCDF file
-        fid = NetCDFFile(self.filename, netcdf_mode_a)
-
-        # Get X, Y from one (any) of the quantities
-        Q = domain.quantities.values()[0]
-        X,Y,_,V = Q.get_vertex_values(xy=True, precision=self.precision)
-
-        # store the connectivity data
-        points = num.concatenate((X[:,num.newaxis],Y[:,num.newaxis]), axis=1)
-        self.writer.store_triangulation(fid,
-                                        points,
-                                        V.astype(num.float32),
-                                        points_georeference=\
-                                            domain.geo_reference)
-
-
-        # Get names of static quantities
-        static_quantities = {}
-        for name in self.writer.static_quantities:
-            Q = domain.quantities[name]
-            A, _ = Q.get_vertex_values(xy=False, 
-                                       precision=self.precision)
-            static_quantities[name] = A
-        
-        # Store static quantities        
-        self.writer.store_static_quantities(fid, **static_quantities)
-                                            
-        fid.close()
-
-    ##
-    # @brief Store time and time dependent quantities 
-    # to the underlying data file.
-    def store_timestep(self):
-        """Store time and time dependent quantities
-        """
-
-        from Scientific.IO.NetCDF import NetCDFFile
-        import types
-        from time import sleep
-        from os import stat
-
-        # Get NetCDF
-        retries = 0
-        file_open = False
-        while not file_open and retries < 10:
-            try:
-                # Open existing file
-                fid = NetCDFFile(self.filename, netcdf_mode_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 step %s again' % self.domain.time
-                log.critical(msg)
-                retries += 1
-                sleep(1)
-            else:
-                file_open = True
-
-        if not file_open:
-            msg = 'File %s could not be opened for append' % self.filename
-            raise DataFileNotOpenError, msg
-
-        # Check to see if the file is already too big:
-        time = fid.variables['time']
-        i = len(time) + 1
-        file_size = stat(self.filename)[6]
-        file_size_increase = file_size / i
-        if file_size + file_size_increase > self.max_size * 2**self.recursion:
-            # In order to get the file name and start time correct,
-            # I change the domain.filename and domain.starttime.
-            # This is the only way to do this without changing
-            # other modules (I think).
-
-            # Write a filename addon that won't break the anuga viewers
-            # (10.sww is bad)
-            filename_ext = '_time_%s' % self.domain.time
-            filename_ext = filename_ext.replace('.', '_')
-
-            # Remember the old filename, then give domain a
-            # name with the extension
-            old_domain_filename = self.domain.get_name()
-            if not self.recursion:
-                self.domain.set_name(old_domain_filename + filename_ext)
-
-            # Temporarily change the domain starttime to the current time
-            old_domain_starttime = self.domain.starttime
-            self.domain.starttime = self.domain.get_time()
-
-            # Build a new data_structure.
-            next_data_structure = SWW_file(self.domain, mode=self.mode,
-                                           max_size=self.max_size,
-                                           recursion=self.recursion+1)
-            if not self.recursion:
-                log.critical('    file_size = %s' % file_size)
-                log.critical('    saving file to %s'
-                             % next_data_structure.filename) 
-
-            # Set up the new data_structure
-            self.domain.writer = next_data_structure
-
-            # Store connectivity and first timestep
-            next_data_structure.store_connectivity()
-            next_data_structure.store_timestep()
-            fid.sync()
-            fid.close()
-
-            # Restore the old starttime and filename
-            self.domain.starttime = old_domain_starttime
-            self.domain.set_name(old_domain_filename)
-        else:
-            self.recursion = False
-            domain = self.domain
-
-            # Get the variables
-            time = fid.variables['time']
-            i = len(time)
-             
-            if 'stage' in self.writer.dynamic_quantities:            
-                # Select only those values for stage, 
-                # xmomentum and ymomentum (if stored) where 
-                # depth exceeds minimum_storable_height
-                #
-                # In this branch it is assumed that elevation
-                # is also available as a quantity            
-            
-                Q = domain.quantities['stage']
-                w, _ = Q.get_vertex_values(xy=False)
-                
-                Q = domain.quantities['elevation']
-                z, _ = Q.get_vertex_values(xy=False)                
-                
-                storable_indices = (w-z >= self.minimum_storable_height)
-            else:
-                # Very unlikely branch
-                storable_indices = None # This means take all
-            
-            
-            # Now store dynamic quantities
-            dynamic_quantities = {}
-            for name in self.writer.dynamic_quantities:
-                netcdf_array = fid.variables[name]
-                
-                Q = domain.quantities[name]
-                A, _ = Q.get_vertex_values(xy=False,
-                                           precision=self.precision)
-
-                if storable_indices is not None:
-                    if name == 'stage':
-                        A = num.choose(storable_indices, (z, A))
-
-                    if name in ['xmomentum', 'ymomentum']:
-                        # Get xmomentum where depth exceeds 
-                        # minimum_storable_height
-                        
-                        # Define a zero vector of same size and type as A
-                        # for use with momenta
-                        null = num.zeros(num.size(A), A.dtype.char)
-                        A = num.choose(storable_indices, (null, A))
-                
-                dynamic_quantities[name] = A
-                
-                                        
-            # Store dynamic quantities
-            self.writer.store_quantities(fid,
-                                         time=self.domain.time,
-                                         sww_precision=self.precision,
-                                         **dynamic_quantities)
-
-
-            # Update extrema if requested
-            domain = self.domain
-            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'][:] = \
-                                        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'][:] = \
-                                        info['max_location']
-                        fid.variables[q + '.max_time'][0] = info['max_time']
-
-            # Flush and close
-            fid.sync()
-            fid.close()
-
-
-##
-# @brief Class to open an sww file so that domain can be populated with quantity values 
-class Read_sww:
-
-    def __init__(self, source):
-        """The source parameter is assumed to be a NetCDF sww file.
-        """
-
-        self.source = source
-
-        self.frame_number = 0
-
-        fin = NetCDFFile(self.source, 'r')
-
-        self.time = num.array(fin.variables['time'], num.float)
-        self.last_frame_number = self.time.shape[0] - 1
-
-        self.frames = num.arange(self.last_frame_number+1)
-
-        fin.close()
-        
-        self.read_mesh()
-
-        self.quantities = {}
-
-        self.read_quantities()
-
-
-    def read_mesh(self):
-        fin = NetCDFFile(self.source, 'r')
-
-        self.vertices = num.array(fin.variables['volumes'], num.int)
-        
-        self.x = x = num.array(fin.variables['x'], num.float)
-        self.y = y = num.array(fin.variables['y'], num.float)
-
-        assert len(self.x) == len(self.y)
-        
-        self.xmin = num.min(x)
-        self.xmax = num.max(x)
-        self.ymin = num.min(y)
-        self.ymax = num.max(y)
-
-
-
-        fin.close()
-        
-    def read_quantities(self, frame_number=0):
-
-        assert frame_number >= 0 and frame_number <= self.last_frame_number
-
-        self.frame_number = frame_number
-
-        M = len(self.x)/3
-        
-        fin = NetCDFFile(self.source, 'r')
-        
-        for q in filter(lambda n:n != 'x' and n != 'y' and n != 'time' and n != 'volumes' and \
-                        '_range' not in n, \
-                        fin.variables.keys()):
-            if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
-                self.quantities[q] = num.ravel(num.array(fin.variables[q], num.float)).reshape(M,3)
-            else: # Time-varying, get the current timestep data
-                self.quantities[q] = num.array(fin.variables[q][self.frame_number], num.float).reshape(M,3)
-        fin.close()
-        return self.quantities
-
-    def get_bounds(self):
-        return [self.xmin, self.xmax, self.ymin, self.ymax]
-
-    def get_last_frame_number(self):
-        return self.last_frame_number
-
-    def get_time(self):
-        return self.time[self.frame_number]
-
-
-##
-# @brief Class for handling checkpoints data
-# @note This is not operational at the moment
-class CPT_file(Data_format):
-    """Interface to native NetCDF format (.cpt) to be 
-    used for checkpointing (one day)
-    """
-
-    ##
-    # @brief Initialize this instantiation.
-    # @param domain ??
-    # @param mode Mode of underlying data file (default WRITE).
-    def __init__(self, domain, mode=netcdf_mode_w):
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        self.precision = netcdf_float #Use full precision
-
-        Data_format.__init__(self, domain, 'sww', mode)
-
-        # NetCDF file definition
-        fid = NetCDFFile(self.filename, mode)
-        if mode[0] == 'w':
-            # Create new file
-            fid.institution = 'Geoscience Australia'
-            fid.description = 'Checkpoint data'
-            #fid.smooth = domain.smooth
-            fid.order = domain.default_order
-
-            # Dimension definitions
-            fid.createDimension('number_of_volumes', self.number_of_volumes)
-            fid.createDimension('number_of_vertices', 3)
-
-            # Store info at all vertices (no smoothing)
-            fid.createDimension('number_of_points', 3*self.number_of_volumes)
-            fid.createDimension('number_of_timesteps', None) #extensible
-
-            # Variable definitions
-
-            # Mesh
-            fid.createVariable('x', self.precision, ('number_of_points',))
-            fid.createVariable('y', self.precision, ('number_of_points',))
-
-
-            fid.createVariable('volumes', netcdf_int, ('number_of_volumes',
-                                                       'number_of_vertices'))
-
-            fid.createVariable('time', self.precision, ('number_of_timesteps',))
-
-            #Allocate space for all quantities
-            for name in domain.quantities.keys():
-                fid.createVariable(name, self.precision,
-                                   ('number_of_timesteps',
-                                    'number_of_points'))
-
-        fid.close()
-
-    ##
-    # @brief Store connectivity data to underlying data file.
-    def store_checkpoint(self):
-        """Write x,y coordinates of triangles.
-        Write connectivity (
-        constituting
-        the bed elevation.
-        """
-
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        domain = self.domain
-
-        #Get NetCDF
-        fid = NetCDFFile(self.filename, netcdf_mode_a)
-
-        # Get the variables
-        x = fid.variables['x']
-        y = fid.variables['y']
-
-        volumes = fid.variables['volumes']
-
-        # Get X, Y and bed elevation Z
-        Q = domain.quantities['elevation']
-        X,Y,Z,V = Q.get_vertex_values(xy=True, precision=self.precision)
-
-        x[:] = X.astype(self.precision)
-        y[:] = Y.astype(self.precision)
-        z[:] = Z.astype(self.precision)
-
-        volumes[:] = V
-
-        fid.close()
-
-    ##
-    # @brief Store time 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
-
-        #Get NetCDF
-        retries = 0
-        file_open = False
-        while not file_open and retries < 10:
-            try:
-                fid = NetCDFFile(self.filename, netcdf_mode_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' \
-                      ' - trying again' % self.filename
-                log.critical(msg)
-                retries += 1
-                sleep(1)
-            else:
-                file_open = True
-
-        if not file_open:
-            msg = 'File %s could not be opened for append' % self.filename
-            raise DataFileNotOpenError, msg
-
-        domain = self.domain
-
-        # Get the variables
-        time = fid.variables['time']
-        stage = fid.variables['stage']
-        i = len(time)
-
-        #Store stage
-        time[i] = self.domain.time
-
-        # Get quantity
-        Q = domain.quantities[name]
-        A,V = Q.get_vertex_values(xy=False, precision=self.precision)
-
-        stage[i,:] = A.astype(self.precision)
-
-        #Flush and close
-        fid.sync()
-        fid.close()
 
 
@@ -924 +184 @@
         #The values are the index positions of file columns.
         self._attribute_dic, self._title_index_dic = \
-            csv2dict(self._file_name, title_check_list=title_check_list)
+            load_csv_as_dict(self._file_name, \
+            title_check_list=title_check_list)
         try:
             #Have code here that handles caps or lower
@@ -1186 +447 @@
     will be excluded 
     
-    See underlying function csv2dict for more details.
-    """
-
-    X, _ = csv2dict(file_name)
+    See underlying function load_csv_as_dict for more details.
+    """
+
+    X, _ = load_csv_as_dict(file_name)
 
     msg = 'Polygon csv file must have 3 or 4 columns'
@@ -1268 +529 @@
 # @brief Convert CSV file to a dictionary of arrays.
 # @param file_name The path to the file to read.
-def csv2array(file_name):
+def load_csv_as_array(file_name):
     """
     Convert CSV files of the form:
@@ -1280 +541 @@
 
 
-    See underlying function csv2dict for more details.
-    """
-
-    X, _ = csv2dict(file_name)
+    See underlying function load_csv_as_dict for more details.
+    """
+
+    X, _ = load_csv_as_dict(file_name)
 
     Y = {}
@@ -1291 +552 @@
     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 dictionary, title as key and column info as value.
-    Also, create a dictionary, title as key and column index as value,
-    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]
-    """
-
-    # FIXME(Ole): Consider dealing with files without headers
-    # FIXME(Ole): Consider a wrapper automatically converting text fields
-    #             to the right type by trying for: int, float, string
-    
-    attribute_dic = {}
-    title_index_dic = {}
-    titles_stripped = [] # List of titles
-
-    reader = csv.reader(file(file_name))
-
-    # Read in and manipulate the title info
-    titles = reader.next()
-    for i, title in enumerate(titles):
-        header = title.strip()
-        titles_stripped.append(header)
-        title_index_dic[header] = i
-    title_count = len(titles_stripped)
-
-    # Check required columns
-    if title_check_list is not None:
-        for title_check in title_check_list:
-            if not title_index_dic.has_key(title_check):
-                msg = 'Reading error. This row is not present %s' % title_check
-                raise IOError, msg
-
-    # Create a dictionary of column values, indexed by column title
-    for line in reader:
-        n = len(line) # Number of entries
-        if n != title_count:
-            msg = 'Entry in file %s had %d columns ' % (file_name, n)
-            msg += 'although there were %d headers' % title_count
-            raise IOError, msg
-        for i, value in enumerate(line):
-            attribute_dic.setdefault(titles_stripped[i], []).append(value)
-
-    return attribute_dic, title_index_dic
 
 
@@ -1390 +595 @@
 
     outfile.close()
-
-
-#########################################################
-#Conversion routines
-########################################################
-
-##
-# @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
-
-    # Get NetCDF
-    FN = create_filename('.', basefilename, 'sww', size)
-    log.critical('Reading from %s' % FN)
-    fid = NetCDFFile(FN, netcdf_mode_r)  #Open existing file for read
-
-    # Get the variables
-    x = fid.variables['x']
-    y = fid.variables['y']
-    z = fid.variables['elevation']
-    time = fid.variables['time']
-    stage = fid.variables['stage']
-
-    M = size  #Number of lines
-    xx = num.zeros((M,3), num.float)
-    yy = num.zeros((M,3), num.float)
-    zz = num.zeros((M,3), num.float)
-
-    for i in range(M):
-        for j in range(3):
-            xx[i,j] = x[i+j*M]
-            yy[i,j] = y[i+j*M]
-            zz[i,j] = z[i+j*M]
-
-    # Write obj for bathymetry
-    FN = create_filename('.', basefilename, 'obj', size)
-    write_obj(FN,xx,yy,zz)
-
-    # Now read all the data with variable information, combine with
-    # x,y info and store as obj
-    for k in range(len(time)):
-        t = time[k]
-        log.critical('Processing timestep %f' % t)
-
-        for i in range(M):
-            for j in range(3):
-                zz[i,j] = stage[k,i+j*M]
-
-        #Write obj for variable data
-        #FN = create_filename(basefilename, 'obj', size, time=t)
-        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
-        write_obj(FN, xx, yy, zz)
-
-
-##
-# @brief 
-# @param basefilename Stem of filename, needs size and extension added.
-def dat2obj(basefilename):
-    """Convert line based data output to obj
-    FIXME: Obsolete?
-    """
-
-    import glob, os
-    from anuga.config import data_dir
-
-    # Get bathymetry and x,y's
-    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
-
-    M = len(lines)  #Number of lines
-    x = num.zeros((M,3), num.float)
-    y = num.zeros((M,3), num.float)
-    z = num.zeros((M,3), num.float)
-
-    for i, line in enumerate(lines):
-        tokens = line.split()
-        values = map(float, tokens)
-
-        for j in range(3):
-            x[i,j] = values[j*3]
-            y[i,j] = values[j*3+1]
-            z[i,j] = values[j*3+2]
-
-    # 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:
-        log.critical('Processing %s' % filename)
-
-        lines = open(data_dir + os.sep + filename, 'r').readlines()
-        assert len(lines) == M
-        root, ext = os.path.splitext(filename)
-
-        # Get time from filename
-        i0 = filename.find('_time=')
-        if i0 == -1:
-            #Skip bathymetry file
-            continue
-
-        i0 += 6  #Position where time starts
-        i1 = filename.find('.dat')
-
-        if i1 > i0:
-            t = float(filename[i0:i1])
-        else:
-            raise DataTimeError, 'Hmmmm'
-
-        for i, line in enumerate(lines):
-            tokens = line.split()
-            values = map(float,tokens)
-
-            for j in range(3):
-                z[i,j] = values[j]
-
-        # Write obj for variable data
-        write_obj(data_dir + os.sep + basefilename + '_time=%.4f' % t, x, y, z)
-
 
 ##
@@ -1573 +653 @@
     outfile.close()
 
-
-##
-# @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,
-            northing_min=None, northing_max=None,
-            use_cache=False, verbose=False,):
-    """Read Digitial Elevation model from the following NetCDF format (.dem)
-
-    Example:
-
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    Convert to NetCDF pts format which is
-
-    points:  (Nx2) float array
-    elevation: N float array
-    """
-
-    kwargs = {'basename_out': basename_out,
-              'easting_min': easting_min,
-              'easting_max': easting_max,
-              'northing_min': northing_min,
-              'northing_max': northing_max,
-              'verbose': verbose}
-
-    if use_cache is True:
-        from caching import cache
-        result = cache(_dem2pts, basename_in, kwargs,
-                       dependencies = [basename_in + '.dem'],
-                       verbose = verbose)
-
-    else:
-        result = apply(_dem2pts, [basename_in], kwargs)
-
-    return result
-
-
-##
-# @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,
-            northing_min=None, northing_max=None):
-    """Read Digitial Elevation model from the following NetCDF format (.dem)
-
-    Internal function. See public function dem2pts for details.
-    """
-
-    # FIXME: Can this be written feasibly using write_pts?
-
-    import os
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    root = basename_in
-
-    # Get NetCDF
-    infile = NetCDFFile(root + '.dem', netcdf_mode_r) 
-
-    if verbose: log.critical('Reading DEM from %s' % (root + '.dem'))
-
-    ncols = infile.ncols[0]
-    nrows = infile.nrows[0]
-    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
-    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
-    cellsize = infile.cellsize[0]
-    NODATA_value = infile.NODATA_value[0]
-    dem_elevation = infile.variables['elevation']
-
-    zone = infile.zone[0]
-    false_easting = infile.false_easting[0]
-    false_northing = infile.false_northing[0]
-
-    # Text strings
-    projection = infile.projection
-    datum = infile.datum
-    units = infile.units
-
-    # Get output file
-    if basename_out == None:
-        ptsname = root + '.pts'
-    else:
-        ptsname = basename_out + '.pts'
-
-    if verbose: log.critical('Store to NetCDF file %s' % ptsname)
-
-    # NetCDF file definition
-    outfile = NetCDFFile(ptsname, netcdf_mode_w)
-
-    # Create new file
-    outfile.institution = 'Geoscience Australia'
-    outfile.description = 'NetCDF pts format for compact and portable ' \
-                          'storage of spatial point data'
-
-    # Assign default values
-    if easting_min is None: easting_min = xllcorner
-    if easting_max is None: easting_max = xllcorner + ncols*cellsize
-    if northing_min is None: northing_min = yllcorner
-    if northing_max is None: northing_max = yllcorner + nrows*cellsize
-
-    # Compute offsets to update georeferencing
-    easting_offset = xllcorner - easting_min
-    northing_offset = yllcorner - northing_min
-
-    # Georeferencing
-    outfile.zone = zone
-    outfile.xllcorner = easting_min # Easting of lower left corner
-    outfile.yllcorner = northing_min # Northing of lower left corner
-    outfile.false_easting = false_easting
-    outfile.false_northing = false_northing
-
-    outfile.projection = projection
-    outfile.datum = datum
-    outfile.units = units
-
-    # Grid info (FIXME: probably not going to be used, but heck)
-    outfile.ncols = ncols
-    outfile.nrows = nrows
-
-    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
-    totalnopoints = nrows*ncols
-
-    # Calculating number of NODATA_values for each row in clipped region
-    # FIXME: use array operations to do faster
-    nn = 0
-    k = 0
-    i1_0 = 0
-    j1_0 = 0
-    thisj = 0
-    thisi = 0
-    for i in range(nrows):
-        y = (nrows-i-1)*cellsize + yllcorner
-        for j in range(ncols):
-            x = j*cellsize + xllcorner
-            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 k == 0:
-                    i1_0 = i
-                    j1_0 = j
-
-                k += 1
-
-    index1 = j1_0
-    index2 = thisj
-
-    # Dimension definitions
-    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
-    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
-
-    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
-    nopoints = clippednopoints-nn
-
-    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
-
-    if verbose:
-        log.critical('There are %d values in the elevation' % totalnopoints)
-        log.critical('There are %d values in the clipped elevation'
-                     % clippednopoints)
-        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
-
-    outfile.createDimension('number_of_points', nopoints)
-    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-
-    # Variable definitions
-    outfile.createVariable('points', netcdf_float, ('number_of_points',
-                                                    'number_of_dimensions'))
-    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
-
-    # Get handles to the variables
-    points = outfile.variables['points']
-    elevation = outfile.variables['elevation']
-
-    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:
-            log.critical('Processing row %d of %d' % (i, nrows))
-
-        lower_index = global_index
-
-        v = dem_elevation_r[i,index1:index2+1]
-        no_NODATA = num.sum(v == NODATA_value)
-        if no_NODATA > 0:
-            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
-        else:
-            newcols = lenv              # ncols_in_bounding_box
-
-        telev = num.zeros(newcols, num.float)
-        tpoints = num.zeros((newcols, 2), num.float)
-
-        local_index = 0
-
-        y = (nrows-i-1)*cellsize + yllcorner
-        #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]
-                telev[local_index] = dem_elevation_r[i, j]
-                global_index += 1
-                local_index += 1
-
-        upper_index = global_index
-
-        if upper_index == lower_index + newcols:
-            points[lower_index:upper_index, :] = tpoints
-            elevation[lower_index:upper_index] = telev
-
-    assert global_index == nopoints, 'index not equal to number of points'
-
-    infile.close()
-    outfile.close()
 
 
@@ -3456 +2293 @@
 
     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:
-
-        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
-
-    E.g
-
-      31/08/04 00:00:00, 1.328223 0 0
-      31/08/04 00:15:00, 1.292912 0 0
-
-    or time (seconds), value0 value1 value2 ...
-
-      0.0, 1.328223 0 0
-      0.1, 1.292912 0 0
-
-    will provide a time dependent function f(t) with three attributes
-
-    filename is assumed to be the rootname with extenisons .txt and .sww
-    """
-
-    import time, calendar
-    from anuga.config import time_format
-    from anuga.utilities.numerical_tools import ensure_numeric
-
-    file_text = filename + '.txt'
-    fid = open(file_text)
-    line = fid.readline()
-    fid.close()
-
-    fields = line.split(',')
-    msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \
-          % file_text
-    assert len(fields) == 2, msg
-
-    if not time_as_seconds:
-        try:
-            starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:
-            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:
-        try:
-            starttime = float(fields[0])
-        except Error:
-            msg = "Bad time format"
-            raise DataTimeError, msg
-
-    # Split values
-    values = []
-    for value in fields[1].split():
-        values.append(float(value))
-
-    q = ensure_numeric(values)
-
-    msg = 'ERROR: File must contain at least one independent value'
-    assert len(q.shape) == 1, msg
-
-    # Read times proper
-    from anuga.config import time_format
-    import time, calendar
-
-    fid = open(file_text)
-    lines = fid.readlines()
-    fid.close()
-
-    N = len(lines)
-    d = len(q)
-
-    T = num.zeros(N, num.float)       # Time
-    Q = num.zeros((N, d), num.float)  # Values
-
-    for i, line in enumerate(lines):
-        fields = line.split(',')
-        if not time_as_seconds:
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-        else:
-             realtime = float(fields[0])
-        T[i] = realtime - starttime
-
-        for j, value in enumerate(fields[1].split()):
-            Q[i, j] = float(value)
-
-    msg = 'File %s must list time as a monotonuosly ' % filename
-    msg += 'increasing sequence'
-    assert num.alltrue(T[1:] - T[:-1] > 0), msg
-
-    #Create NetCDF file
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
-
-    fid.institution = 'Geoscience Australia'
-    fid.description = 'Time series'
-
-    #Reference point
-    #Start time in seconds since the epoch (midnight 1/1/1970)
-    #FIXME: Use Georef
-    fid.starttime = starttime
-
-    # dimension definitions
-    #fid.createDimension('number_of_volumes', self.number_of_volumes)
-    #fid.createDimension('number_of_vertices', 3)
-
-    fid.createDimension('number_of_timesteps', len(T))
-
-    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
-
-    fid.variables['time'][:] = T
-
-    for i in range(Q.shape[1]):
-        try:
-            name = quantity_names[i]
-        except:
-            name = 'Attribute%d' % i
-
-        fid.createVariable(name, netcdf_float, ('number_of_timesteps',))
-        fid.variables[name][:] = Q[:,i]
-
-    fid.close()
 
 
@@ -4054 +2760 @@
 
 ##
-# @brief 
-# @param filename
-# @param verbose
-def tsh2sww(filename, verbose=False):
-    """
-    to check if a tsh/msh file 'looks' good.
-    """
-
-    if verbose == True: log.critical('Creating domain from %s' % filename)
-
-    domain = pmesh_to_domain_instance(filename, Domain)
-
-    if verbose == True: log.critical("Number of triangles = %s" % len(domain))
-
-    domain.smooth = True
-    domain.format = 'sww'   #Native netcdf visualisation format
-    file_path, filename = path.split(filename)
-    filename, ext = path.splitext(filename)
-    domain.set_name(filename)
-    domain.reduction = mean
-
-    if verbose == True: log.critical("file_path = %s" % file_path)
-
-    if file_path == "":
-        file_path = "."
-    domain.set_datadir(file_path)
-
-    if verbose == True:
-        log.critical("Output written to %s%s%s.%s"
-                     % (domain.get_datadir(), sep, domain.get_name(),
-                        domain.format))
-
-    sww = SWW_file(domain)
-    sww.store_connectivity()
-    sww.store_timestep()
-
-
-##
-# @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,
-                  ucur_dir,
-                  vcur_dir,
-                  sww_file,
-                  minlat=None, maxlat=None,
-                  minlon=None, maxlon=None,
-                  zscale=1,
-                  mean_stage=0,
-                  fail_on_NaN=True,
-                  elevation_NaN_filler=0,
-                  bath_prefix='ba',
-                  elevation_prefix='el',
-                  verbose=False):
-    """
-    Produce an sww boundary file, from esri ascii data from CSIRO.
-
-    Also convert latitude and longitude to UTM. All coordinates are
-    assumed to be given in the GDA94 datum.
-
-    assume:
-    All files are in esri ascii format
-
-    4 types of information
-    bathymetry
-    elevation
-    u velocity
-    v velocity
-
-    Assumptions
-    The metadata of all the files is the same
-    Each type is in a seperate directory
-    One bath file with extention .000
-    The time period is less than 24hrs and uniform.
-    """
-
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    from anuga.coordinate_transforms.redfearn import redfearn
-
-    precision = netcdf_float # So if we want to change the precision its done here
-
-    # go in to the bath dir and load the only file,
-    bath_files = os.listdir(bath_dir)
-    bath_file = bath_files[0]
-    bath_dir_file =  bath_dir + os.sep + bath_file
-    bath_metadata, bath_grid =  _read_asc(bath_dir_file)
-
-    #Use the date.time of the bath file as a basis for
-    #the start time for other files
-    base_start = bath_file[-12:]
-
-    #go into the elevation dir and load the 000 file
-    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
-                         + base_start
-
-    elevation_files = os.listdir(elevation_dir)
-    ucur_files = os.listdir(ucur_dir)
-    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
-    assert elevation_files[0] == 'el' + base_start
-
-    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'] \
-                  for x in range(number_of_longitudes)]
-    latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \
-                 for y in range(number_of_latitudes)]
-
-     # reverse order of lat, so the 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)
-
-    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
-    latitudes = latitudes[kmin:kmax]
-    longitudes = longitudes[lmin:lmax]
-    number_of_latitudes = len(latitudes)
-    number_of_longitudes = len(longitudes)
-    number_of_times = len(os.listdir(elevation_dir))
-    number_of_points = number_of_latitudes * number_of_longitudes
-    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
-
-    #Work out the times
-    if len(elevation_files) > 1:
-        # Assume: The time period is less than 24hrs.
-        time_period = (int(elevation_files[1][-3:]) \
-                       - int(elevation_files[0][-3:])) * 60*60
-        times = [x*time_period for x in range(len(elevation_files))]
-    else:
-        times = [0.0]
-
-    if verbose:
-        log.critical('------------------------------------------------')
-        log.critical('Statistics:')
-        log.critical('  Extent (lat/lon):')
-        log.critical('    lat in [%f, %f], len(lat) == %d'
-                     % (min(latitudes), max(latitudes), len(latitudes)))
-        log.critical('    lon in [%f, %f], len(lon) == %d'
-                     % (min(longitudes), max(longitudes), len(longitudes)))
-        log.critical('    t in [%f, %f], len(t) == %d'
-                     % (min(times), max(times), len(times)))
-
-    ######### WRITE THE SWW FILE #############
-
-    # NetCDF file definition
-    outfile = NetCDFFile(sww_file, netcdf_mode_w)
-
-    #Create new file
-    outfile.institution = 'Geoscience Australia'
-    outfile.description = 'Converted from XXX'
-
-    #For sww compatibility
-    outfile.smoothing = 'Yes'
-    outfile.order = 1
-
-    #Start time in seconds since the epoch (midnight 1/1/1970)
-    outfile.starttime = starttime = times[0]
-
-    # dimension definitions
-    outfile.createDimension('number_of_volumes', number_of_volumes)
-    outfile.createDimension('number_of_vertices', 3)
-    outfile.createDimension('number_of_points', number_of_points)
-    outfile.createDimension('number_of_timesteps', number_of_times)
-
-    # variable definitions
-    outfile.createVariable('x', precision, ('number_of_points',))
-    outfile.createVariable('y', precision, ('number_of_points',))
-    outfile.createVariable('elevation', precision, ('number_of_points',))
-
-    #FIXME: Backwards compatibility
-    #outfile.createVariable('z', precision, ('number_of_points',))
-    #################################
-
-    outfile.createVariable('volumes', netcdf_int, ('number_of_volumes',
-                                                   'number_of_vertices'))
-
-    outfile.createVariable('time', precision, ('number_of_timesteps',))
-
-    outfile.createVariable('stage', precision, ('number_of_timesteps',
-                                                'number_of_points'))
-
-    outfile.createVariable('xmomentum', precision, ('number_of_timesteps',
-                                                    'number_of_points'))
-
-    outfile.createVariable('ymomentum', precision, ('number_of_timesteps',
-                                                    'number_of_points'))
-
-    #Store
-    from anuga.coordinate_transforms.redfearn import redfearn
-
-    x = num.zeros(number_of_points, num.float)  #Easting
-    y = num.zeros(number_of_points, num.float)  #Northing
-
-    if verbose: log.critical('Making triangular grid')
-
-    #Get zone of 1st point.
-    refzone, _, _ = redfearn(latitudes[0], longitudes[0])
-
-    vertices = {}
-    i = 0
-    for k, lat in enumerate(latitudes):
-        for l, lon in enumerate(longitudes):
-            vertices[l,k] = i
-
-            zone, easting, northing = redfearn(lat, lon)
-
-            #msg = 'Zone boundary crossed at longitude =', lon
-            #assert zone == refzone, msg
-            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
-            x[i] = easting
-            y[i] = northing
-            i += 1
-
-    #Construct 2 triangles per 'rectangular' element
-    volumes = []
-    for l in range(number_of_longitudes-1):    #X direction
-        for k in range(number_of_latitudes-1): #Y direction
-            v1 = vertices[l,k+1]
-            v2 = vertices[l,k]
-            v3 = vertices[l+1,k+1]
-            v4 = vertices[l+1,k]
-
-            #Note, this is different to the ferrit2sww code
-            #since the order of the lats is reversed.
-            volumes.append([v1,v3,v2]) #Upper element
-            volumes.append([v4,v2,v3]) #Lower element
-
-    volumes = num.array(volumes, num.int)      #array default#
-
-    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., (max(x)+min(y))/2.)
-
-    if verbose:
-        log.critical('------------------------------------------------')
-        log.critical('More Statistics:')
-        log.critical('  Extent (/lon):')
-        log.critical('    x in [%f, %f], len(lat) == %d'
-                     % (min(x), max(x), len(x)))
-        log.critical('    y in [%f, %f], len(lon) == %d'
-                     % (min(y), max(y), len(y)))
-        log.critical('geo_ref: ', geo_ref)
-
-    z = num.resize(bath_grid,outfile.variables['elevation'][:].shape)
-    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
-    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
-    # 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(num.int32) # On Opteron 64
-
-    stage = outfile.variables['stage']
-    xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum']
-
-    outfile.variables['time'][:] = times   #Store time relative
-
-    if verbose: log.critical('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])
-
-        #cut matrix to desired size
-        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
-        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
-        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
-
-        # handle missing values
-        missing = (elevation_grid == elevation_meta['NODATA_value'])
-        if num.sometrue (missing):
-            if fail_on_NaN:
-                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: log.critical('  Doing %d of %d'
-                                                          % (j, n))
-
-        i = 0
-        for k in range(number_of_latitudes):      #Y direction
-            for l in range(number_of_longitudes): #X direction
-                w = zscale*elevation_grid[k,l] + mean_stage
-                stage[j,i] = w
-                h = w - z[i]
-                xmomentum[j,i] = u_momentum_grid[k,l]*h
-                ymomentum[j,i] = v_momentum_grid[k,l]*h
-                i += 1
-
-    outfile.close()
-
-
-##
 # @brief Return max&min indexes (for slicing) of area specified.
 # @param latitudes_ref ??
@@ -5951 +4328 @@
 
 
-# @brief A class to write an SWW file.
-class Write_sww:
-    from anuga.shallow_water.shallow_water_domain import Domain
-
-    RANGE = '_range'
-    EXTREMA = ':extrema'
-
-    ##
-    # brief Instantiate the SWW writer class.
-    def __init__(self, static_quantities, dynamic_quantities):
-        """Initialise Write_sww with two list af quantity names: 
-        
-        static_quantities (e.g. elevation or friction): 
-            Stored once at the beginning of the simulation in a 1D array
-            of length number_of_points   
-        dynamic_quantities (e.g stage):
-            Stored every timestep in a 2D array with 
-            dimensions number_of_points X number_of_timesteps        
-        
-        """
-        self.static_quantities = static_quantities   
-        self.dynamic_quantities = dynamic_quantities
-
-
-    ##
-    # @brief Store a header in the SWW file.
-    # @param outfile Open handle to the file that will be written.
-    # @param times A list of time slices *or* a start time.
-    # @param number_of_volumes The number of triangles.
-    # @param number_of_points The number of points.
-    # @param description The internal file description string.
-    # @param smoothing True if smoothing is to be used.
-    # @param order 
-    # @param sww_precision Data type of the quantity written (netcdf constant)
-    # @param verbose True if this function is to be verbose.
-    # @note If 'times' is a list, the info will be made relative.
-    def store_header(self,
-                     outfile,
-                     times,
-                     number_of_volumes,
-                     number_of_points,
-                     description='Generated by ANUGA',
-                     smoothing=True,
-                     order=1,
-                     sww_precision=netcdf_float32,
-                     verbose=False):
-        """Write an SWW file header.
-
-        outfile - the open file that will be written
-        times - A list of the time slice times OR a start time
-        Note, if a list is given the info will be made relative.
-        number_of_volumes - the number of triangles
-        """
-
-        outfile.institution = 'Geoscience Australia'
-        outfile.description = description
-
-        # For sww compatibility
-        if smoothing is True:
-            # Smoothing to be depreciated
-            outfile.smoothing = 'Yes'
-            outfile.vertices_are_stored_uniquely = 'False'
-        else:
-            # Smoothing to be depreciated
-            outfile.smoothing = 'No'
-            outfile.vertices_are_stored_uniquely = 'True'
-        outfile.order = order
-
-        try:
-            revision_number = get_revision_number()
-        except:
-            revision_number = None
-        # Allow None to be stored as a string
-        outfile.revision_number = str(revision_number)
-
-        # This is being used to seperate one number from a list.
-        # what it is actually doing is sorting lists from numeric arrays.
-        if isinstance(times, (list, num.ndarray)):
-            number_of_times = len(times)
-            times = ensure_numeric(times)
-            if number_of_times == 0:
-                starttime = 0
-            else:
-                starttime = times[0]
-                times = times - starttime  #Store relative times
-        else:
-            number_of_times = 0
-            starttime = times
-
-
-        outfile.starttime = starttime
-
-        # dimension definitions
-        outfile.createDimension('number_of_volumes', number_of_volumes)
-        outfile.createDimension('number_of_vertices', 3)
-        outfile.createDimension('numbers_in_range', 2)
-
-        if smoothing is True:
-            outfile.createDimension('number_of_points', number_of_points)
-            # FIXME(Ole): This will cause sww files for parallel domains to
-            # have ghost nodes stored (but not used by triangles).
-            # To clean this up, we have to change get_vertex_values and
-            # friends in quantity.py (but I can't be bothered right now)
-        else:
-            outfile.createDimension('number_of_points', 3*number_of_volumes)
-
-        outfile.createDimension('number_of_timesteps', number_of_times)
-
-        # variable definitions
-        outfile.createVariable('x', sww_precision, ('number_of_points',))
-        outfile.createVariable('y', sww_precision, ('number_of_points',))
-
-        outfile.createVariable('volumes', netcdf_int, ('number_of_volumes',
-                                                       'number_of_vertices'))
-
-        # Doing sww_precision instead of Float gives cast errors.
-        outfile.createVariable('time', netcdf_float,
-                               ('number_of_timesteps',))
-
-                               
-        for q in self.static_quantities:
-            
-            outfile.createVariable(q, sww_precision,
-                                   ('number_of_points',))
-            
-            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
-                                   ('numbers_in_range',))
-                                   
-            # Initialise ranges with small and large sentinels.
-            # If this was in pure Python we could have used None sensibly
-            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
-            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
-
-        #if 'elevation' in self.static_quantities:    
-        #    # FIXME: Backwards compat - get rid of z once old view has retired
-        #    outfile.createVariable('z', sww_precision,
-        #                           ('number_of_points',))
-                               
-        for q in self.dynamic_quantities:
-            outfile.createVariable(q, sww_precision, ('number_of_timesteps',
-                                                      'number_of_points'))
-            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
-                                   ('numbers_in_range',))
-
-            # Initialise ranges with small and large sentinels.
-            # If this was in pure Python we could have used None sensibly
-            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
-            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
-
-        if isinstance(times, (list, num.ndarray)):
-            outfile.variables['time'][:] = times    # Store time relative
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('Statistics:')
-            log.critical('    t in [%f, %f], len(t) == %d'
-                         % (num.min(times), num.max(times), len(times.flat)))
-
-    ##
-    # @brief Store triangulation data in the underlying file.
-    # @param outfile Open handle to underlying file.
-    # @param points_utm List or array of points in UTM.
-    # @param volumes 
-    # @param zone 
-    # @param new_origin georeference that the points can be set to.
-    # @param points_georeference The georeference of the points_utm.
-    # @param verbose True if this function is to be verbose.
-    def store_triangulation(self,
-                            outfile,
-                            points_utm,
-                            volumes,
-                            zone=None, 
-                            new_origin=None,
-                            points_georeference=None, 
-                            verbose=False):
-        """
-        new_origin - qa georeference that the points can be set to. (Maybe
-        do this before calling this function.)
-
-        points_utm - currently a list or array of the points in UTM.
-        points_georeference - the georeference of the points_utm
-
-        How about passing new_origin and current_origin.
-        If you get both, do a convertion from the old to the new.
-
-        If you only get new_origin, the points are absolute,
-        convert to relative
-
-        if you only get the current_origin the points are relative, store
-        as relative.
-
-        if you get no georefs create a new georef based on the minimums of
-        points_utm.  (Another option would be to default to absolute)
-
-        Yes, and this is done in another part of the code.
-        Probably geospatial.
-
-        If you don't supply either geo_refs, then supply a zone. If not
-        the default zone will be used.
-
-        precon:
-            header has been called.
-        """
-
-        number_of_points = len(points_utm)
-        volumes = num.array(volumes)
-        points_utm = num.array(points_utm)
-
-        # Given the two geo_refs and the points, do the stuff
-        # described in the method header
-        # if this is needed else where, pull out as a function
-        points_georeference = ensure_geo_reference(points_georeference)
-        new_origin = ensure_geo_reference(new_origin)
-        if new_origin is None and points_georeference is not None:
-            points = points_utm
-            geo_ref = points_georeference
-        else:
-            if new_origin is None:
-                new_origin = Geo_reference(zone, min(points_utm[:,0]),
-                                                 min(points_utm[:,1]))
-            points = new_origin.change_points_geo_ref(points_utm,
-                                                      points_georeference)
-            geo_ref = new_origin
-
-        # At this stage I need a georef and points
-        # the points are relative to the georef
-        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.)
-
-        x =  points[:,0]
-        y =  points[:,1]
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('More Statistics:')
-            log.critical('  Extent (/lon):')
-            log.critical('    x in [%f, %f], len(lat) == %d'
-                         % (min(x), max(x), len(x)))
-            log.critical('    y in [%f, %f], len(lon) == %d'
-                         % (min(y), max(y), len(y)))
-            #log.critical('    z in [%f, %f], len(z) == %d'
-            #             % (min(elevation), max(elevation), len(elevation)))
-            log.critical('geo_ref: %s' % str(geo_ref))
-            log.critical('------------------------------------------------')
-
-        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
-        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
-        outfile.variables['volumes'][:] = volumes.astype(num.int32) #On Opteron 64
-
-
-
-    # @brief Write the static quantity data to the underlying file.
-    # @param outfile Handle to open underlying file.
-    # @param sww_precision Format of quantity data to write (default Float32).
-    # @param verbose True if this function is to be verbose.
-    # @param **quant
-    def store_static_quantities(self, 
-                                outfile, 
-                                sww_precision=num.float32,
-                                verbose=False, 
-                                **quant):
-        """
-        Write the static quantity info.
-
-        **quant is extra keyword arguments passed in. These must be
-          the numpy arrays to be stored in the sww file at each timestep.
-
-        The argument sww_precision allows for storing as either 
-        * single precision (default): num.float32
-        * double precision: num.float64 or num.float 
-
-        Precondition:
-            store_triangulation and
-            store_header have been called.
-        """
-
-        # The dictionary quant must contain numpy arrays for each name.
-        # These will typically be quantities from Domain such as friction 
-        #
-        # Arrays not listed in static_quantitiues will be ignored, silently.
-        #
-        # This method will also write the ranges for each quantity, 
-        # e.g. stage_range, xmomentum_range and ymomentum_range
-        for q in self.static_quantities:
-            if not quant.has_key(q):
-                msg = 'Values for quantity %s was not specified in ' % q
-                msg += 'store_quantities so they cannot be stored.'
-                raise NewQuantity, msg
-            else:
-                q_values = ensure_numeric(quant[q])
-                
-                x = q_values.astype(sww_precision)
-                outfile.variables[q][:] = x
-        
-                # This populates the _range values
-                outfile.variables[q + Write_sww.RANGE][0] = num.min(x)
-                outfile.variables[q + Write_sww.RANGE][1] = num.max(x)
-                    
-        # FIXME: Hack for backwards compatibility with old viewer
-        #if 'elevation' in self.static_quantities:
-        #    outfile.variables['z'][:] = outfile.variables['elevation'][:]
-
-                    
-                    
-        
-        
-    ##
-    # @brief Write the quantity data to the underlying file.
-    # @param outfile Handle to open underlying file.
-    # @param sww_precision Format of quantity data to write (default Float32).
-    # @param slice_index
-    # @param time
-    # @param verbose True if this function is to be verbose.
-    # @param **quant
-    def store_quantities(self, 
-                         outfile, 
-                         sww_precision=num.float32,
-                         slice_index=None,
-                         time=None,
-                         verbose=False, 
-                         **quant):
-        """
-        Write the quantity info at each timestep.
-
-        **quant is extra keyword arguments passed in. These must be
-          the numpy arrays to be stored in the sww file at each timestep.
-
-        if the time array is already been built, use the slice_index
-        to specify the index.
-
-        Otherwise, use time to increase the time dimension
-
-        Maybe make this general, but the viewer assumes these quantities,
-        so maybe we don't want it general - unless the viewer is general
-        
-        The argument sww_precision allows for storing as either 
-        * single precision (default): num.float32
-        * double precision: num.float64 or num.float 
-
-        Precondition:
-            store_triangulation and
-            store_header have been called.
-        """
-
-        if time is not None:
-            file_time = outfile.variables['time']
-            slice_index = len(file_time)
-            file_time[slice_index] = time
-        else:
-            slice_index = int(slice_index) # Has to be cast in case it was numpy.int    
-
-        # Write the named dynamic quantities
-        # The dictionary quant must contain numpy arrays for each name.
-        # These will typically be the conserved quantities from Domain 
-        # (Typically stage,  xmomentum, ymomentum).
-        #
-        # Arrays not listed in dynamic_quantitiues will be ignored, silently.
-        #
-        # This method will also write the ranges for each quantity, 
-        # e.g. stage_range, xmomentum_range and ymomentum_range
-        for q in self.dynamic_quantities:
-            if not quant.has_key(q):
-                msg = 'Values for quantity %s was not specified in ' % q
-                msg += 'store_quantities so they cannot be stored.'
-                raise NewQuantity, msg
-            else:
-                q_values = ensure_numeric(quant[q])
-                
-                x = q_values.astype(sww_precision)
-                outfile.variables[q][slice_index] = x
-                    
-        
-                # This updates the _range values
-                q_range = outfile.variables[q + Write_sww.RANGE][:]
-                q_values_min = num.min(q_values)
-                if q_values_min < q_range[0]:
-                    outfile.variables[q + Write_sww.RANGE][0] = q_values_min
-                q_values_max = num.max(q_values)
-                if q_values_max > q_range[1]:
-                    outfile.variables[q + Write_sww.RANGE][1] = q_values_max
-
-    ##
-    # @brief Print the quantities in the underlying file.
-    # @param outfile UNUSED.
-    def verbose_quantities(self, outfile):
-        log.critical('------------------------------------------------')
-        log.critical('More Statistics:')
-        for q in self.dynamic_quantities:
-            log.critical('  %s in [%f, %f]'
-                         % (q, outfile.variables[q+Write_sww.RANGE][0],
-                            outfile.variables[q+Write_sww.RANGE][1]))
-        log.critical('------------------------------------------------')
-
-
 ##
 # @brief Obsolete?
@@ -7652 +5633 @@
 
 ##
-# @brief Find all SWW files in a directory with given stem name.
-# @param look_in_dir The directory to look in.
-# @param base_name The file stem name.
-# @param verbose True if this function is to be verbose.
-# @return A list of found filename strings.
-# @note Will accept 'base_name' with or without '.sww' extension.
-# @note If no files found, raises IOError exception.
-def get_all_swwfiles(look_in_dir='', base_name='', verbose=False):
-    '''
-    Finds all the sww files in a "look_in_dir" which contains a "base_name".
-    will accept base_name with or without the extension ".sww"
-
-    Returns: a list of strings
-
-    Usage:     iterate_over = get_all_swwfiles(dir, name)
-    then
-               for swwfile in iterate_over:
-                   do stuff
-
-    Check "export_grids" and "get_maximum_inundation_data" for examples
-    '''
-
-    # plus tests the extension
-    name, extension = os.path.splitext(base_name)
-
-    if extension != '' and extension != '.sww':
-        msg = 'file %s%s must be a NetCDF sww file!' % (base_name, extension)
-        raise IOError, msg
-
-    if look_in_dir == "":
-        look_in_dir = "."                                   # Unix compatibility
-
-    dir_ls = os.listdir(look_in_dir)
-    iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == '.sww']
-    if len(iterate_over) == 0:
-        msg = 'No files of the base name %s' % name
-        raise IOError, msg
-
-    if verbose: log.critical('iterate over %s' % iterate_over)
-
-    return iterate_over
-
-
-##
-# @brief Find all files in a directory that contain a string and have extension.
-# @param look_in_dir Path to the directory to look in.
-# @param base_name Stem filename of the file(s) of interest.
-# @param extension Extension of the files to look for.
-# @param verbose True if this function is to be verbose.
-# @return A list of found filename strings.
-# @note If no files found, raises IOError exception.
-def get_all_files_with_extension(look_in_dir='',
-                                 base_name='',
-                                 extension='.sww',
-                                 verbose=False):
-    '''Find all files in a directory with given stem name.
-    Finds all the sww files in a "look_in_dir" which contains a "base_name".
-
-    Returns: a list of strings
-
-    Usage:     iterate_over = get_all_swwfiles(dir, name)
-    then
-               for swwfile in iterate_over:
-                   do stuff
-
-    Check "export_grids" and "get_maximum_inundation_data" for examples
-    '''
-
-    # plus tests the extension
-    name, ext = os.path.splitext(base_name)
-
-    if ext != '' and ext != extension:
-        msg = 'base_name %s must be a file with %s extension!' \
-              % (base_name, extension)
-        raise IOError, msg
-
-    if look_in_dir == "":
-        look_in_dir = "."                               # Unix compatibility
-
-    dir_ls = os.listdir(look_in_dir)
-    iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == extension]
-
-    if len(iterate_over) == 0:
-        msg = 'No files of the base name %s in %s' % (name, look_in_dir)
-        raise IOError, msg
-
-    if verbose: log.critical('iterate over %s' % iterate_over)
-
-    return iterate_over
-
-
-##
-# @brief Find all files in a directory that contain a given string.
-# @param look_in_dir Path to the directory to look in.
-# @param base_name String that files must contain.
-# @param verbose True if this function is to be verbose.
-def get_all_directories_with_name(look_in_dir='', base_name='', verbose=False):
-    '''
-    Finds all the directories in a "look_in_dir" which contains a "base_name".
-
-    Returns: a list of strings
-
-    Usage:     iterate_over = get_all_directories_with_name(dir, name)
-    then:      for swwfile in iterate_over:
-                   do stuff
-
-    Check "export_grids" and "get_maximum_inundation_data" for examples
-    '''
-
-    if look_in_dir == "":
-        look_in_dir = "."                                  # Unix compatibility
-
-    dir_ls = os.listdir(look_in_dir)
-    iterate_over = [x for x in dir_ls if base_name in x]
-
-    if len(iterate_over) == 0:
-        msg = 'No files of the base name %s' % base_name
-        raise IOError, msg
-
-    if verbose: log.critical('iterate over %s' % iterate_over)
-
-    return iterate_over
-
-
-##
 # @brief Convert points to a polygon (?)
 # @param points_file The points file.

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -664 +664 @@
         """
 
-        from anuga.shallow_water.data_manager import SWW_file
+        from anuga.shallow_water.sww_file import SWW_file
         
         # Initialise writer

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -18 +18 @@
 from anuga.shallow_water import *
 from anuga.shallow_water.data_manager import *
+from anuga.shallow_water.sww_file import SWW_file
+from anuga.shallow_water.file_conversion import tsh2sww, \
+                        asc_csiro2sww, pmesh_to_domain_instance, \
+                        dem2pts
 from anuga.config import epsilon, g
 from anuga.utilities.anuga_exceptions import ANUGAError
@@ -24 +28 @@
 from anuga.abstract_2d_finite_volumes.util import file_function
 from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.utilities.file_utils import del_dir, load_csv_as_dict
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.config import netcdf_float
@@ -4444 +4449 @@
 
         import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
 
         # the memory optimised least squares
@@ -4469 +4473 @@
 
         import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
 
         from data_manager import _read_asc
@@ -7710 +7713 @@
             x = fid.variables['x'][:]+fid.xllcorner    # x-coordinates of vertices
             y = fid.variables['y'][:]+fid.yllcorner    # y-coordinates of vertices
-            elevation = fid.variables['elevation'][:]
+            elevation = fid.variables['elevation'][:]
             time=fid.variables['time'][:]+fid.starttime
 
@@ -7736 +7739 @@
                 count = 0
 
-                # read in urs test data for stage, e and n velocity
+                # read in urs test data for stage, e and n velocity
                 urs_file_name_z = 'z_'+str(source_number)+'_'+str(j)+'.csv'
-                dict = csv2array(os.path.join(testdir, urs_file_name_z))
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
                 urs_stage = dict['urs_stage']
                 urs_file_name_e = 'e_'+str(source_number)+'_'+str(j)+'.csv'
-                dict = csv2array(os.path.join(testdir, urs_file_name_e))
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
                 urs_e = dict['urs_e']
                 urs_file_name_n = 'n_'+str(source_number)+'_'+str(j)+'.csv'
-                dict = csv2array(os.path.join(testdir, urs_file_name_n))
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
                 urs_n = dict['urs_n']
 
-                # find start and end time for stage             
+                # find start and end time for stage             
                 for i in range(len(urs_stage)):
                     if urs_stage[i] == 0.0:
@@ -7774 +7777 @@
                 assert num.allclose(index_start_urs_z,start_times_z[j]/delta_t), msg
 
-                msg = 'e velocity start time from urs file is not the same as the '
+                msg = 'e velocity start time from urs file is not the same as the '
                 msg += 'header file for source %i and station %i' %(source_number,j)
-                assert num.allclose(index_start_urs_e,start_times_e[j]/delta_t), msg
-
-                msg = 'n velocity start time from urs file is not the same as the '
+                assert num.allclose(index_start_urs_e,start_times_e[j]/delta_t), msg
+
+                msg = 'n velocity start time from urs file is not the same as the '
                 msg += 'header file for source %i and station %i' %(source_number,j)
-                assert num.allclose(index_start_urs_n,start_times_n[j]/delta_t), msg
-                
-                # get index for start and end time for sts quantities
+                assert num.allclose(index_start_urs_n,start_times_n[j]/delta_t), msg
+                
+                # get index for start and end time for sts quantities
                 index_start_stage = 0
                 index_end_stage = 0
@@ -7810 +7813 @@
                                     sts_stage[index_start_stage:index_end_stage], 
                                     rtol=1.0e-6, atol=1.0e-5 ), msg                                
-                                
+                                
                 # check that urs e velocity and sts xmomentum are the same
-                msg = 'urs e velocity is not equivalent to sts x momentum for for source %i and station %i' %(source_number,j)
+                msg = 'urs e velocity is not equivalent to sts x momentum for for source %i and station %i' %(source_number,j)
                 assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]),
                                 sts_xmom[index_start_x:index_end_x], 
                                 rtol=1.0e-5, atol=1.0e-4 ), msg
-                
+                
                 # check that urs n velocity and sts ymomentum are the same
                 #print 'urs n velocity', urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j])
@@ -7864 +7867 @@
         x = fid.variables['x'][:]+fid.xllcorner   # x-coordinates of vertices
         y = fid.variables['y'][:]+fid.yllcorner   # y-coordinates of vertices
-        elevation = fid.variables['elevation'][:]
+        elevation = fid.variables['elevation'][:]
         time=fid.variables['time'][:]+fid.starttime
         
@@ -7892 +7895 @@
         assert num.allclose(sts_starttime, starttime), msg
     
-        #stations = [1,2,3]
+        #stations = [1,2,3]
         #for j in stations: 
         for j in range(len(x)):
             index_start_urs_z = 0
             index_end_urs_z = 0
-            index_start_urs_e = 0
+            index_start_urs_e = 0
             index_end_urs_e = 0
-            index_start_urs_n = 0
+            index_start_urs_n = 0
             index_end_urs_n = 0
             count = 0
 
-            # read in urs test data for stage, e and n velocity
+            # read in urs test data for stage, e and n velocity
             urs_file_name_z = 'z_combined_'+str(j)+'.csv'
-            dict = csv2array(os.path.join(testdir, urs_file_name_z))
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
             urs_stage = dict['urs_stage']
             urs_file_name_e = 'e_combined_'+str(j)+'.csv'
-            dict = csv2array(os.path.join(testdir, urs_file_name_e))
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
             urs_e = dict['urs_e']
             urs_file_name_n = 'n_combined_'+str(j)+'.csv'
-            dict = csv2array(os.path.join(testdir, urs_file_name_n))
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
             urs_n = dict['urs_n']
 
-            # find start and end time for stage         
+            # find start and end time for stage         
             for i in range(len(urs_stage)):
                 if urs_stage[i] == 0.0:
@@ -7937 +7940 @@
             assert num.allclose(index_start_urs_z,start_times_z/delta_t), msg
 
-            msg = 'e velocity start time from urs file is not the same as the '
+            msg = 'e velocity start time from urs file is not the same as the '
             msg += 'header file at station %i' %(j)
-            assert num.allclose(index_start_urs_e,start_times_e/delta_t), msg
-
-            msg = 'n velocity start time from urs file is not the same as the '
+            assert num.allclose(index_start_urs_e,start_times_e/delta_t), msg
+
+            msg = 'n velocity start time from urs file is not the same as the '
             msg += 'header file at station %i' %(j)
-            assert num.allclose(index_start_urs_n,start_times_n/delta_t), msg
-                
-            # get index for start and end time for sts quantities
+            assert num.allclose(index_start_urs_n,start_times_n/delta_t), msg
+                
+            # get index for start and end time for sts quantities
             index_start_stage = 0
             index_end_stage = 0
@@ -7963 +7966 @@
                     index_end_stage = i
                 
-            index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
-
-            index_start_x = index_start_stage
-            index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e])
-            sts_xmom = quantities['ymomentum'][:,j]
-
-            index_start_y = index_start_stage
-            index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n])
-            sts_ymom = quantities['ymomentum'][:,j]
+            index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
+
+            index_start_x = index_start_stage
+            index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e])
+            sts_xmom = quantities['ymomentum'][:,j]
+
+            index_start_y = index_start_stage
+            index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n])
+            sts_ymom = quantities['ymomentum'][:,j]
 
             # check that urs stage and sts stage are the same
             msg = 'urs stage is not equal to sts stage for station %i' %j
-            #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z]
+            #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z]
             #print 'sts stage', sts_stage[index_start_stage:index_end_stage]
             #print 'diff', max(urs_stage[index_start_urs_z:index_end_urs_z]-sts_stage[index_start_stage:index_end_stage])
@@ -7982 +7985 @@
                             sts_stage[index_start_stage:index_end_stage], 
                                 rtol=1.0e-5, atol=1.0e-4 ), msg                                
-                                
+                                
             # check that urs e velocity and sts xmomentum are the same          
             msg = 'urs e velocity is not equivalent to sts xmomentum for station %i' %j
@@ -7988 +7991 @@
                             sts_xmom[index_start_x:index_end_x], 
                             rtol=1.0e-5, atol=1.0e-4 ), msg
-                
+                
             # check that urs n velocity and sts ymomentum are the same                            
             msg = 'urs n velocity is not equivalent to sts ymomentum for station %i' %j
@@ -7995 +7998 @@
                             rtol=1.0e-5, atol=1.0e-4 ), msg
 
-        fid.close()
+        fid.close()
         
         os.remove(sts_name_out+'.sts')

anuga_core/source/anuga/shallow_water/test_forcing.py

@@ -10 +10 @@
 import numpy as num
 
-### Helper functions
-###
 
 def scalar_func_list(t, x, y):
@@ -21 +19 @@
     return [17.7]
 
-# Variable windfield implemented using functions
+
 def speed(t, x, y):
-    """Large speeds halfway between center and edges
+    """
+    Variable windfield implemented using functions    
+    Large speeds halfway between center and edges
 
     Low speeds at center and edges
@@ -248 +248 @@
 
         # Convert ASCII file to NetCDF (Which is what we really like!)
-        from data_manager import timefile2netcdf
+        from file_conversion import timefile2netcdf
 
         timefile2netcdf(filename)
@@ -340 +340 @@
 
         # Convert ASCII file to NetCDF (Which is what we really like!)
-        from data_manager import timefile2netcdf
+        from file_conversion import timefile2netcdf
 
         timefile2netcdf(filename, time_as_seconds=True)

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -6070 +6070 @@
         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
         from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
+        from anuga.shallow_water.boundaries import Reflective_boundary
         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
         from anuga.shallow_water.shallow_water_domain import Rainfall

anuga_core/source/anuga/shallow_water/tsh2sww.py

@@ -13 +13 @@
 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
 import time, os 
-from anuga.pyvolution.data_manager import SWW_file
+from anuga.pyvolution.sww_file import SWW_file
 from anuga.utilities.numerical_tools import mean
 import anuga.utilities.log as log