Index: /anuga_core/source/anuga/file_conversion/sww2dem.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/sww2dem.py	(revision 7744)
+++ /anuga_core/source/anuga/file_conversion/sww2dem.py	(revision 7744)
@@ -0,0 +1,415 @@
+# external modules
+import numpy as num
+
+# ANUGA modules
+from anuga.abstract_2d_finite_volumes.util import remove_lone_verts     
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.utilities.system_tools import get_vars_in_expression
+import anuga.utilities.log as log
+
+
+######
+# formula mappings
+######
+
+quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
+                    'depth':'stage-elevation',
+                    'speed': \
+ '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
+
+
+
+# Default block size for sww2dem()
+DEFAULT_BLOCK_SIZE = 10000
+
+def sww2dem(basename_in, basename_out=None,
+            quantity=None, # defaults to elevation
+            reduction=None,
+            cellsize=10,
+            number_of_decimal_places=None,
+            NODATA_value=-9999,
+            easting_min=None,
+            easting_max=None,
+            northing_min=None,
+            northing_max=None,
+            verbose=False,
+            origin=None,
+            datum='WGS84',
+            format='ers',
+            block_size=None):
+    """Read SWW file and convert to Digitial Elevation model format
+    (.asc or .ers)
+
+    Example (ASC):
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    The number of decimal places can be specified by the user to save
+    on disk space requirements by specifying in the call to sww2dem.
+
+    Also write accompanying file with same basename_in but extension .prj
+    used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
+
+
+    The parameter quantity must be the name of an existing quantity or
+    an expression involving existing quantities. The default is
+    'elevation'. Quantity is not a list of quantities.
+
+    if timestep (an index) is given, output quantity at that timestep
+
+    if reduction is given and its an index, output quantity at that timestep. If reduction is given
+    and is a built in function, use that to reduce quantity over all timesteps.
+
+    datum
+
+    format can be either 'asc' or 'ers'
+    block_size - sets the number of slices along the non-time axis to
+                 process in one block.
+    """
+
+    import sys
+    import types
+
+    from anuga.geometry.polygon import inside_polygon, outside_polygon
+    from anuga.abstract_2d_finite_volumes.util import \
+         apply_expression_to_dictionary
+
+    msg = 'Format must be either asc or ers'
+    assert format.lower() in ['asc', 'ers'], msg
+
+    false_easting = 500000
+    false_northing = 10000000
+
+    if quantity is None:
+        quantity = 'elevation'
+    
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' % quantity
+
+    if quantity_formula.has_key(quantity):
+        quantity = quantity_formula[quantity]
+
+    if number_of_decimal_places is None:
+        number_of_decimal_places = 3
+
+    if block_size is None:
+        block_size = DEFAULT_BLOCK_SIZE
+
+    # Read SWW file
+    swwfile = basename_in + '.sww'
+    demfile = basename_out + '.' + format
+
+    # Read sww file
+    if verbose:
+        log.critical('Reading from %s' % swwfile)
+        log.critical('Output directory is %s' % basename_out)
+
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    #Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+    if type(reduction) is not types.BuiltinFunctionType:
+        times = fid.variables['time'][reduction]
+    else:
+        times = fid.variables['time'][:]
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+
+    if origin is None:
+        # Get geo_reference
+        # sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() # Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+    # FIXME: Refactor using code from Interpolation_function.statistics
+    # (in interpolate.py)
+    # Something like print swwstats(swwname)
+    if verbose:
+        log.critical('------------------------------------------------')
+        log.critical('Statistics of SWW file:')
+        log.critical('  Name: %s' % swwfile)
+        log.critical('  Reference:')
+        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
+        if type(reduction) is not types.BuiltinFunctionType:
+            log.critical('    Time: %f' % times)
+        else:
+            log.critical('    Start time: %f' % fid.starttime[0])
+        log.critical('  Extent:')
+        log.critical('    x [m] in [%f, %f], len(x) == %d'
+                     %(num.min(x), num.max(x), len(x.flat)))
+        log.critical('    y [m] in [%f, %f], len(y) == %d'
+                     % (num.min(y), num.max(y), len(y.flat)))
+        if type(reduction) is not types.BuiltinFunctionType:
+            log.critical('    t [s] = %f, len(t) == %d' % (times, 1))
+        else:
+            log.critical('    t [s] in [%f, %f], len(t) == %d'
+                         % (min(times), max(times), len(times)))
+        log.critical('  Quantities [SI units]:')
+        
+        # Comment out for reduced memory consumption
+        for name in ['stage', 'xmomentum', 'ymomentum']:
+            q = fid.variables[name][:].flatten()
+            if type(reduction) is not types.BuiltinFunctionType:
+                q = q[reduction*len(x):(reduction+1)*len(x)]
+            if verbose: log.critical('    %s in [%f, %f]'
+                                     % (name, min(q), max(q)))
+        for name in ['elevation']:
+            q = fid.variables[name][:].flatten()
+            if verbose: log.critical('    %s in [%f, %f]'
+                                     % (name, min(q), max(q)))
+
+    # Get the variables in the supplied expression.
+    # This may throw a SyntaxError exception.
+    var_list = get_vars_in_expression(quantity)
+
+    # Check that we have the required variables in the SWW file.
+    missing_vars = []
+    for name in var_list:
+        try:
+            _ = fid.variables[name]
+        except:
+            missing_vars.append(name)
+    if missing_vars:
+        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
+               % (quantity, swwfile))
+        raise Exception, msg
+
+    # Create result array and start filling, block by block.
+    result = num.zeros(number_of_points, num.float)
+
+    for start_slice in xrange(0, number_of_points, block_size):
+        # Limit slice size to array end if at last block
+        end_slice = min(start_slice + block_size, number_of_points)
+        
+        # Get slices of all required variables
+        q_dict = {}
+        for name in var_list:
+            # check if variable has time axis
+            if len(fid.variables[name].shape) == 2:
+                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
+            else:       # no time axis
+                q_dict[name] = fid.variables[name][start_slice:end_slice]
+
+        # Evaluate expression with quantities found in SWW file
+        res = apply_expression_to_dictionary(quantity, q_dict)
+
+        if len(res.shape) == 2:
+            new_res = num.zeros(res.shape[1], num.float)
+            for k in xrange(res.shape[1]):
+                if type(reduction) is not types.BuiltinFunctionType:
+                    new_res[k] = res[reduction,k]
+                else:
+                    new_res[k] = reduction(res[:,k])
+            res = new_res
+
+        result[start_slice:end_slice] = res
+                                    
+    # Post condition: Now q has dimension: number_of_points
+    assert len(result.shape) == 1
+    assert result.shape[0] == number_of_points
+
+    if verbose:
+        log.critical('Processed values for %s are in [%f, %f]'
+                     % (quantity, min(result), max(result)))
+
+    # Create grid and update xll/yll corner and x,y
+    # Relative extent
+    if easting_min is None:
+        xmin = min(x)
+    else:
+        xmin = easting_min - xllcorner
+
+    if easting_max is None:
+        xmax = max(x)
+    else:
+        xmax = easting_max - xllcorner
+
+    if northing_min is None:
+        ymin = min(y)
+    else:
+        ymin = northing_min - yllcorner
+
+    if northing_max is None:
+        ymax = max(y)
+    else:
+        ymax = northing_max - yllcorner
+
+    msg = 'xmax must be greater than or equal to xmin.\n'
+    msg += 'I got xmin = %f, xmax = %f' %(xmin, xmax)
+    assert xmax >= xmin, msg
+
+    msg = 'ymax must be greater than or equal to xmin.\n'
+    msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
+    assert ymax >= ymin, msg
+
+    if verbose: log.critical('Creating grid')
+    ncols = int((xmax-xmin)/cellsize) + 1
+    nrows = int((ymax-ymin)/cellsize) + 1
+
+    # New absolute reference and coordinates
+    newxllcorner = xmin + xllcorner
+    newyllcorner = ymin + yllcorner
+
+    x = x + xllcorner - newxllcorner
+    y = y + yllcorner - newyllcorner
+
+    vertex_points = num.concatenate ((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
+    assert len(vertex_points.shape) == 2
+
+    grid_points = num.zeros ((ncols*nrows, 2), num.float)
+
+    for i in xrange(nrows):
+        if format.lower() == 'asc':
+            yg = i * cellsize
+        else:
+            # this will flip the order of the y values for ers
+            yg = (nrows-i) * cellsize
+
+        for j in xrange(ncols):
+            xg = j * cellsize
+            k = i*ncols + j
+
+            grid_points[k, 0] = xg
+            grid_points[k, 1] = yg
+
+    # Interpolate
+    from anuga.fit_interpolate.interpolate import Interpolate
+
+    # Remove loners from vertex_points, volumes here
+    vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
+    # export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
+    interp = Interpolate(vertex_points, volumes, verbose = verbose)
+
+    # Interpolate using quantity values
+    if verbose: log.critical('Interpolating')
+    grid_values = interp.interpolate(result, grid_points).flatten()
+
+    if verbose:
+        log.critical('Interpolated values are in [%f, %f]'
+                     % (num.min(grid_values), num.max(grid_values)))
+
+    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
+    P = interp.mesh.get_boundary_polygon()
+    outside_indices = outside_polygon(grid_points, P, closed=True)
+
+    for i in outside_indices:
+        grid_values[i] = NODATA_value
+
+    if format.lower() == 'ers':
+        # setup ERS header information
+        grid_values = num.reshape(grid_values, (nrows, ncols))
+        header = {}
+        header['datum'] = '"' + datum + '"'
+        # FIXME The use of hardwired UTM and zone number needs to be made optional
+        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
+        header['projection'] = '"UTM-' + str(zone) + '"'
+        header['coordinatetype'] = 'EN'
+        if header['coordinatetype'] == 'LL':
+            header['longitude'] = str(newxllcorner)
+            header['latitude'] = str(newyllcorner)
+        elif header['coordinatetype'] == 'EN':
+            header['eastings'] = str(newxllcorner)
+            header['northings'] = str(newyllcorner)
+        header['nullcellvalue'] = str(NODATA_value)
+        header['xdimension'] = str(cellsize)
+        header['ydimension'] = str(cellsize)
+        header['value'] = '"' + quantity + '"'
+        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
+
+        #Write
+        if verbose: log.critical('Writing %s' % demfile)
+
+        import ermapper_grids
+
+        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
+
+        fid.close()
+    else:
+        #Write to Ascii format
+        #Write prj file
+        prjfile = basename_out + '.prj'
+
+        if verbose: log.critical('Writing %s' % prjfile)
+        prjid = open(prjfile, 'w')
+        prjid.write('Projection    %s\n' %'UTM')
+        prjid.write('Zone          %d\n' %zone)
+        prjid.write('Datum         %s\n' %datum)
+        prjid.write('Zunits        NO\n')
+        prjid.write('Units         METERS\n')
+        prjid.write('Spheroid      %s\n' %datum)
+        prjid.write('Xshift        %d\n' %false_easting)
+        prjid.write('Yshift        %d\n' %false_northing)
+        prjid.write('Parameters\n')
+        prjid.close()
+
+        if verbose: log.critical('Writing %s' % demfile)
+
+        ascid = open(demfile, 'w')
+
+        ascid.write('ncols         %d\n' %ncols)
+        ascid.write('nrows         %d\n' %nrows)
+        ascid.write('xllcorner     %d\n' %newxllcorner)
+        ascid.write('yllcorner     %d\n' %newyllcorner)
+        ascid.write('cellsize      %f\n' %cellsize)
+        ascid.write('NODATA_value  %d\n' %NODATA_value)
+
+        #Get bounding polygon from mesh
+        #P = interp.mesh.get_boundary_polygon()
+        #inside_indices = inside_polygon(grid_points, P)
+
+        # change printoptions so that a long string of zeros in not
+        # summarized as [0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0]
+        #printoptions = num.get_printoptions()
+        #num.set_printoptions(threshold=sys.maxint)
+
+        format = '%.'+'%g' % number_of_decimal_places +'e'
+        for i in range(nrows):
+            if verbose and i % ((nrows+10)/10) == 0:
+                log.critical('Doing row %d of %d' % (i, nrows))
+
+            base_index = (nrows-i-1)*ncols
+
+            slice = grid_values[base_index:base_index+ncols]
+
+            num.savetxt(ascid, slice.reshape(1,ncols), format, ' ' )
+            
+        
+        #Close
+        ascid.close()
+        fid.close()
+
+        return basename_out
Index: /anuga_core/source/anuga/file_conversion/test_file_conversion.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7743)
+++ /anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7744)
@@ -2241,4 +2241,102 @@
 
 
+    def test_read_asc(self):
+        """Test conversion from dem in ascii format to native NetCDF format
+        """
+
+        import time, os
+
+        from file_conversion import _read_asc
+        #Write test asc file
+        filename = tempfile.mktemp(".000")
+        fid = open(filename, 'w')
+        fid.write("""ncols         7
+nrows         4
+xllcorner     2000.5
+yllcorner     3000.5
+cellsize      25
+NODATA_value  -9999
+    97.921    99.285   125.588   180.830   258.645   342.872   415.836
+   473.157   514.391   553.893   607.120   678.125   777.283   883.038
+   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
+   502.645   516.230   504.739   450.604   388.500   338.097   514.980
+""")
+        fid.close()
+        bath_metadata, grid = _read_asc(filename, verbose=self.verbose)
+        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
+        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
+        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
+        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
+        self.failUnless(grid[0][0]  == 97.921,  'Failed')
+        self.failUnless(grid[3][6]  == 514.980,  'Failed')
+
+        os.remove(filename)
+
+
+    #### TESTS FOR URS 2 SWW  ###     
+    
+    def create_mux(self, points_num=None):
+        # write all the mux stuff.
+        time_step_count = 3
+        time_step = 0.5
+        
+        longitudes = [150.66667, 150.83334, 151., 151.16667]
+        latitudes = [-34.5, -34.33333, -34.16667, -34]
+
+        if points_num == None:
+            points_num = len(longitudes) * len(latitudes)
+
+        lonlatdeps = []
+        quantities = ['HA','UA','VA']
+        mux_names = [WAVEHEIGHT_MUX_LABEL,
+                     EAST_VELOCITY_LABEL,
+                     NORTH_VELOCITY_LABEL]
+        quantities_init = [[],[],[]]
+        # urs binary is latitude fastest
+        for i,lon in enumerate(longitudes):
+            for j,lat in enumerate(latitudes):
+                _ , e, n = redfearn(lat, lon)
+                lonlatdeps.append([lon, lat, n])
+                quantities_init[0].append(e) # HA
+                quantities_init[1].append(n ) # UA
+                quantities_init[2].append(e) # VA
+        #print "lonlatdeps",lonlatdeps
+
+        file_handle, base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(base_name)
+        
+        files = []        
+        for i,q in enumerate(quantities): 
+            quantities_init[i] = ensure_numeric(quantities_init[i])
+            #print "HA_init", HA_init
+            q_time = num.zeros((time_step_count, points_num), num.float)
+            for time in range(time_step_count):
+                q_time[time,:] = quantities_init[i] #* time * 4
+            
+            #Write C files
+            columns = 3 # long, lat , depth
+            file = base_name + mux_names[i]
+            #print "base_name file",file 
+            f = open(file, 'wb')
+            files.append(file)
+            f.write(pack('i',points_num))
+            f.write(pack('i',time_step_count))
+            f.write(pack('f',time_step))
+
+            #write lat/long info
+            for lonlatdep in lonlatdeps:
+                for float in lonlatdep:
+                    f.write(pack('f',float))
+                    
+            # Write quantity info
+            for time in  range(time_step_count):
+                for point_i in range(points_num):
+                    f.write(pack('f',q_time[time,point_i]))
+                    #print " mux_names[i]", mux_names[i] 
+                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
+            f.close()
+        return base_name, files
+
 #-------------------------------------------------------------
 
Index: /anuga_core/source/anuga/file_conversion/test_sww2dem.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/test_sww2dem.py	(revision 7744)
+++ /anuga_core/source/anuga/file_conversion/test_sww2dem.py	(revision 7744)
@@ -0,0 +1,1401 @@
+import unittest, os
+import numpy as num
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+from anuga.shallow_water.shallow_water_domain import Domain
+from anuga.config import netcdf_mode_r
+
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     DEFAULT_ZONE
+
+from anuga.shallow_water.sww_file import SWW_file
+     
+# boundary functions
+from anuga.shallow_water.boundaries import Reflective_boundary, \
+            Field_boundary, Transmissive_momentum_set_stage_boundary, \
+            Transmissive_stage_zero_momentum_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary, Dirichlet_boundary, \
+            Time_boundary, File_boundary, AWI_boundary
+
+# local modules
+from sww2dem import sww2dem
+
+class Test_Sww2Dem(unittest.TestCase):
+    def setUp(self):
+        import time
+        from mesh_factory import rectangular
+        
+        # Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        # Set some field values
+        domain.set_quantity('elevation', lambda x,y: -x)
+        domain.set_quantity('friction', 0.03)
+
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        ######################
+        #Initial condition - with jumps
+        bed = domain.quantities['elevation'].vertex_values
+        stage = num.zeros(bed.shape, num.float)
+
+        h = 0.3
+        for i in range(stage.shape[0]):
+            if i % 2 == 0:
+                stage[i,:] = bed[i,:] + h
+            else:
+                stage[i,:] = bed[i,:]
+
+        domain.set_quantity('stage', stage)
+
+
+        domain.distribute_to_vertices_and_edges()               
+        self.domain = domain
+
+        C = domain.get_vertex_coordinates()
+        self.X = C[:,0:6:2].copy()
+        self.Y = C[:,1:6:2].copy()
+
+        self.F = bed
+        self.verbose = False
+
+
+    def tearDown(self):
+        pass
+
+
+    def test_sww2dem_asc_elevation_depth(self):
+        """test_sww2dem_asc_elevation_depth
+        
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        
+        Also check geo_reference is correct
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        # Setup
+        self.domain.set_name('datatest')
+
+        prjfile = self.domain.get_name() + '_elevation.prj'
+        ascfile = self.domain.get_name() + '_elevation.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+        
+        # Check georeferencig: zone, xllcorner and yllcorner
+        assert fid.zone[0] == 56
+        assert fid.xllcorner[0] == 308500
+        assert fid.yllcorner[0] == 6189000
+                
+
+        fid.close()
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                verbose = self.verbose,
+                format = 'asc')
+
+        #Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                quantity = 'depth',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                verbose = self.verbose,
+                format = 'asc')
+        
+        #Check asc file
+        ascfile = self.domain.get_name() + '_depth.asc'
+        prjfile = self.domain.get_name() + '_depth.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+    def test_sww2dem_larger(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView. Here:
+
+        ncols         11
+        nrows         11
+        xllcorner     308500
+        yllcorner     6189000
+        cellsize      10.000000
+        NODATA_value  -9999
+        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
+         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
+         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
+         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
+         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
+         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
+         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
+         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
+         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
+         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
+           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
+
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Create basic mesh (100m x 100m)
+        points, vertices, boundary = rectangular(2, 2, 100, 100)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        domain.set_name('datatest')
+
+        prjfile = domain.get_name() + '_elevation.prj'
+        ascfile = domain.get_name() + '_elevation.asc'
+        swwfile = domain.get_name() + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        #
+        domain.set_quantity('elevation', lambda x,y: -x-y)
+        domain.set_quantity('stage', 0)
+
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        #
+        sww = SWW_file(domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        
+        domain.tight_slope_limiters = 1
+        domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 10  #10m grid
+
+
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                verbose = self.verbose,
+                format = 'asc',
+                block_size=2)
+
+
+        #Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '11'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '11'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
+        for i, line in enumerate(lines[6:]):
+            for j, value in enumerate( line.split() ):
+                assert num.allclose(float(value), -(10-i+j)*cellsize,
+                                    atol=1.0e-12, rtol=1.0e-12)
+
+                # Note: Equality can be obtained in this case,
+                # but it is better to use allclose.
+                #assert float(value) == -(10-i+j)*cellsize
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+    def test_sww2dem_larger_zero(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. Arcview. This example has rows with a
+        large number of zeros
+
+        ncols         2001
+        nrows         2
+        xllcorner     308500
+        yllcorner     6189000
+        cellsize      1.000000
+        NODATA_value  -9999
+        0.0 ....
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+
+        from mesh_factory import rectangular_cross
+
+        #Create basic mesh (100m x 100m)
+        points, vertices, boundary = rectangular_cross(2000, 1, 2000.0, 1.0)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 1
+
+        domain.set_name('datatest')
+
+        prjfile = domain.get_name() + '_elevation.prj'
+        ascfile = domain.get_name() + '_elevation.asc'
+        swwfile = domain.get_name() + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        #
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 0)
+
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        #
+        sww = SWW_file(domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        
+        domain.tight_slope_limiters = 1
+        domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 1.0  #0.1 grid
+
+
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                verbose = self.verbose,
+                format = 'asc',
+                block_size=2)
+
+
+        #Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '2001'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '2'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
+        for i, line in enumerate(lines[6:]):
+            for j, value in enumerate( line.split() ):
+                assert num.allclose(float(value), 0.0,
+                                    atol=1.0e-12, rtol=1.0e-12)
+
+                # Note: Equality can be obtained in this case,
+                # but it is better to use allclose.
+                #assert float(value) == -(10-i+j)*cellsize
+
+
+        fid.close()
+
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+
+    def test_sww2dem_boundingbox(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView.
+        This will test that mesh can be restricted by bounding box
+
+        Original extent is 100m x 100m:
+
+        Eastings:   308500 -  308600
+        Northings: 6189000 - 6189100
+
+        Bounding box changes this to the 50m x 50m square defined by
+
+        Eastings:   308530 -  308570
+        Northings: 6189050 - 6189100
+
+        The cropped values should be
+
+         -130 -140 -150 -160 -170
+         -120 -130 -140 -150 -160
+         -110 -120 -130 -140 -150
+         -100 -110 -120 -130 -140
+          -90 -100 -110 -120 -130
+          -80  -90 -100 -110 -120
+
+        and the new lower reference point should be
+        Eastings:   308530
+        Northings: 6189050
+
+        Original dataset is the same as in test_sww2dem_larger()
+
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+
+        from mesh_factory import rectangular
+
+        #Create basic mesh (100m x 100m)
+        points, vertices, boundary = rectangular(2, 2, 100, 100)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        domain.set_name('datatest')
+
+        prjfile = domain.get_name() + '_elevation.prj'
+        ascfile = domain.get_name() + '_elevation.asc'
+        swwfile = domain.get_name() + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        #
+        domain.set_quantity('elevation', lambda x,y: -x-y)
+        domain.set_quantity('stage', 0)
+
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        #
+        sww = SWW_file(domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        #domain.tight_slope_limiters = 1
+        domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 10  #10m grid
+
+
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        # Export to ascii/prj files
+        sww2dem(domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                easting_min = 308530,
+                easting_max = 308570,
+                northing_min = 6189050,
+                northing_max = 6189100,
+                verbose = self.verbose,
+                format = 'asc')
+
+        fid.close()
+
+
+        # Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '6'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308530)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189050)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values
+        for i, line in enumerate(lines[6:]):
+            for j, value in enumerate( line.split() ):
+                #assert float(value) == -(10-i+j)*cellsize
+                assert float(value) == -(10-i+j+3)*cellsize
+
+
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+    def test_sww2dem_asc_stage_reduction(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This tests the reduction of quantity stage using min
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        prjfile = self.domain.get_name() + '_stage.prj'
+        ascfile = self.domain.get_name() + '_stage.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        #self.domain.tight_slope_limiters = 1
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                quantity = 'stage',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                reduction = min,
+                format = 'asc',
+                verbose=self.verbose)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+
+        #Check grid values (where applicable)
+        for j in range(5):
+            if j%2 == 0:
+                L = lines[6+j].strip().split()
+                jj = 4-j
+                for i in range(5):
+                    if i%2 == 0:
+                        index = jj/2 + i/2*3
+                        val0 = stage[0,index]
+                        val1 = stage[1,index]
+
+                        #print i, j, index, ':', L[i], val0, val1
+                        assert num.allclose(float(L[i]), min(val0, val1))
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_sww2dem_asc_stage_time(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This tests the reduction of quantity stage using min
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        prjfile = self.domain.get_name() + '_stage.prj'
+        ascfile = self.domain.get_name() + '_stage.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        #self.domain.tight_slope_limiters = 1
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                quantity = 'stage',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                reduction = 1,
+                format = 'asc',
+                verbose=self.verbose)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values (where applicable)
+        for j in range(5):
+            if j%2 == 0:
+                L = lines[6+j].strip().split()
+                jj = 4-j
+                for i in range(5):
+                    if i%2 == 0:
+                        index = jj/2 + i/2*3
+                        
+                        val = stage[1,index]
+                   
+                        assert num.allclose(float(L[i]), val)
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+    def test_sww2dem_asc_derived_quantity(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This tests the use of derived quantities
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        prjfile = self.domain.get_name() + '_depth.prj'
+        ascfile = self.domain.get_name() + '_depth.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 0.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        #self.domain.tight_slope_limiters = 1
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                basename_out = 'datatest_depth',
+                quantity = 'stage - elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                reduction = min,
+                format = 'asc',
+                verbose = self.verbose)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+
+        #Check grid values (where applicable)
+        for j in range(5):
+            if j%2 == 0:
+                L = lines[6+j].strip().split()
+                jj = 4-j
+                for i in range(5):
+                    if i%2 == 0:
+                        index = jj/2 + i/2*3
+                        val0 = stage[0,index] - z[index]
+                        val1 = stage[1,index] - z[index]
+
+                        #print i, j, index, ':', L[i], val0, val1
+                        assert num.allclose(float(L[i]), min(val0, val1))
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+
+
+    def test_sww2dem_asc_missing_points(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This test includes the writing of missing values
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup mesh not coinciding with rectangle.
+        #This will cause missing values to occur in gridded data
+
+
+        points = [                        [1.0, 1.0],
+                              [0.5, 0.5], [1.0, 0.5],
+                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
+
+        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
+
+        #Create shallow water domain
+        domain = Domain(points, vertices)
+        domain.default_order=2
+
+
+        #Set some field values
+        domain.set_quantity('elevation', lambda x,y: -x-y)
+        domain.set_quantity('friction', 0.03)
+
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'exterior': B} )
+
+
+        ######################
+        #Initial condition - with jumps
+
+        bed = domain.quantities['elevation'].vertex_values
+        stage = num.zeros(bed.shape, num.float)
+
+        h = 0.3
+        for i in range(stage.shape[0]):
+            if i % 2 == 0:
+                stage[i,:] = bed[i,:] + h
+            else:
+                stage[i,:] = bed[i,:]
+
+        domain.set_quantity('stage', stage)
+        domain.distribute_to_vertices_and_edges()
+
+        domain.set_name('datatest')
+
+        prjfile = domain.get_name() + '_elevation.prj'
+        ascfile = domain.get_name() + '_elevation.asc'
+        swwfile = domain.get_name() + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+
+        domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(swwfile, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
+
+        #Export to ascii/prj files
+        sww2dem(domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                verbose = self.verbose,
+                format = 'asc')
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert num.allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            assert len(L) == 5
+            y = (4-j) * cellsize
+
+            for i in range(5):
+                #print i
+                if i+j >= 4:
+                    assert num.allclose(float(L[i]), -i*cellsize - y)
+                else:
+                    #Missing values
+                    assert num.allclose(float(L[i]), -9999)
+
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+    def test_sww2ers_simple(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+
+        NODATA_value = 1758323
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        headerfile = self.domain.get_name() + '.ers'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+
+        #self.domain.tight_slope_limiters = 1
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ers files
+        sww2dem(self.domain.get_name(),
+                quantity = 'elevation',
+                cellsize = cellsize,
+                number_of_decimal_places = 9,
+                NODATA_value = NODATA_value,
+                verbose = self.verbose,
+                format = 'ers')
+
+        #Check header data
+        from ermapper_grids import read_ermapper_header, read_ermapper_data
+
+        header = read_ermapper_header(self.domain.get_name() + '_elevation.ers')
+        #print header
+        assert header['projection'].lower() == '"utm-56"'
+        assert header['datum'].lower() == '"wgs84"'
+        assert header['units'].lower() == '"meters"'
+        assert header['value'].lower() == '"elevation"'
+        assert header['xdimension'] == '0.25'
+        assert header['ydimension'] == '0.25'
+        assert float(header['eastings']) == 308500.0   #xllcorner
+        assert float(header['northings']) == 6189000.0 #yllcorner
+        assert int(header['nroflines']) == 5
+        assert int(header['nrofcellsperline']) == 5
+        assert int(header['nullcellvalue']) == NODATA_value
+        #FIXME - there is more in the header
+
+
+        #Check grid data
+        grid = read_ermapper_data(self.domain.get_name() + '_elevation')
+
+        #FIXME (Ole): Why is this the desired reference grid for -x-y?
+        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
+                    -1,    -1.25, -1.5,  -1.75, -2.0,
+                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
+                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
+                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
+
+
+        #print grid
+        assert num.allclose(grid, ref_grid)
+
+        fid.close()
+
+        #Cleanup
+        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
+        #Done (Ole) - it was because sww2ers didn't close it's sww file
+        os.remove(sww.filename)
+        os.remove(self.domain.get_name() + '_elevation')
+        os.remove(self.domain.get_name() + '_elevation.ers')
+
+#################################################################################
+
+if __name__ == "__main__":
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
+    suite = unittest.makeSuite(Test_Sww2Dem, 'test')
+    runner = unittest.TextTestRunner(verbosity=1)
+    runner.run(suite)
Index: /anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7743)
+++ /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7744)
@@ -88,4 +88,5 @@
 from anuga.load_mesh.loadASCII import export_mesh_file
 from anuga.geometry.polygon import intersection
+from anuga.file_conversion.sww2dem import sww2dem
 
 from anuga.utilities.system_tools import get_vars_in_expression
@@ -99,20 +100,4 @@
                 DataFileNotOpenError, DataTimeError, DataDomainError, \
                 NewQuantity
-
-
-# Default block size for sww2dem()
-DEFAULT_BLOCK_SIZE = 10000
-
-
-######
-# formula mappings
-######
-
-quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
-                    'depth':'stage-elevation',
-                    'speed': \
- '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
-
-
 
 
@@ -279,30 +264,4 @@
 
             
-            
-##
-# @brief Convert CSV file to a dictionary of arrays.
-# @param file_name The path to the file to read.
-def load_csv_as_array(file_name):
-    """
-    Convert CSV files of the form:
-
-    time, discharge, velocity
-    0.0,  1.2,       0.0
-    0.1,  3.2,       1.1
-    ...
-
-    to a dictionary of numeric arrays.
-
-
-    See underlying function load_csv_as_dict for more details.
-    """
-
-    X, _ = load_csv_as_dict(file_name)
-
-    Y = {}
-    for key in X.keys():
-        Y[key] = num.array([float(x) for x in X[key]])
-
-    return Y
 
 
@@ -764,481 +723,4 @@
                                               verbose=verbose)
 
-
-##
-# @brief Convert SWW file to DEM file (.asc or .ers).
-# @param basename_in 
-# @param basename_out 
-# @param quantity 
-# @param timestep 
-# @param reduction 
-# @param cellsize  
-# @param number_of_decimal_places 
-# @param NODATA_value 
-# @param easting_min 
-# @param easting_max 
-# @param northing_min 
-# @param northing_max 
-# @param verbose 
-# @param origin 
-# @param datum 
-# @param format 
-# @return 
-def sww2dem(basename_in, basename_out=None,
-            quantity=None, # defaults to elevation
-            reduction=None,
-            cellsize=10,
-            number_of_decimal_places=None,
-            NODATA_value=-9999,
-            easting_min=None,
-            easting_max=None,
-            northing_min=None,
-            northing_max=None,
-            verbose=False,
-            origin=None,
-            datum='WGS84',
-            format='ers',
-            block_size=None):
-    """Read SWW file and convert to Digitial Elevation model format
-    (.asc or .ers)
-
-    Example (ASC):
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    The number of decimal places can be specified by the user to save
-    on disk space requirements by specifying in the call to sww2dem.
-
-    Also write accompanying file with same basename_in but extension .prj
-    used to fix the UTM zone, datum, false northings and eastings.
-
-    The prj format is assumed to be as
-
-    Projection    UTM
-    Zone          56
-    Datum         WGS84
-    Zunits        NO
-    Units         METERS
-    Spheroid      WGS84
-    Xshift        0.0000000000
-    Yshift        10000000.0000000000
-    Parameters
-
-
-    The parameter quantity must be the name of an existing quantity or
-    an expression involving existing quantities. The default is
-    'elevation'. Quantity is not a list of quantities.
-
-    if timestep (an index) is given, output quantity at that timestep
-
-    if reduction is given and its an index, output quantity at that timestep. If reduction is given
-    and is a built in function, use that to reduce quantity over all timesteps.
-
-    datum
-
-    format can be either 'asc' or 'ers'
-    block_size - sets the number of slices along the non-time axis to
-                 process in one block.
-    """
-
-    import sys
-    import types
-
-    from anuga.geometry.polygon import inside_polygon, outside_polygon
-    from anuga.abstract_2d_finite_volumes.util import \
-         apply_expression_to_dictionary
-
-    msg = 'Format must be either asc or ers'
-    assert format.lower() in ['asc', 'ers'], msg
-
-    false_easting = 500000
-    false_northing = 10000000
-
-    if quantity is None:
-        quantity = 'elevation'
-    
-    if reduction is None:
-        reduction = max
-
-    if basename_out is None:
-        basename_out = basename_in + '_%s' % quantity
-
-    if quantity_formula.has_key(quantity):
-        quantity = quantity_formula[quantity]
-
-    if number_of_decimal_places is None:
-        number_of_decimal_places = 3
-
-    if block_size is None:
-        block_size = DEFAULT_BLOCK_SIZE
-
-    # Read SWW file
-    swwfile = basename_in + '.sww'
-    demfile = basename_out + '.' + format
-
-    # Read sww file
-    if verbose:
-        log.critical('Reading from %s' % swwfile)
-        log.critical('Output directory is %s' % basename_out)
-
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(swwfile)
-
-    #Get extent and reference
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    volumes = fid.variables['volumes'][:]
-    if type(reduction) is not types.BuiltinFunctionType:
-        times = fid.variables['time'][reduction]
-    else:
-        times = fid.variables['time'][:]
-
-    number_of_timesteps = fid.dimensions['number_of_timesteps']
-    number_of_points = fid.dimensions['number_of_points']
-
-    if origin is None:
-        # Get geo_reference
-        # sww files don't have to have a geo_ref
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference() # Default georef object
-
-        xllcorner = geo_reference.get_xllcorner()
-        yllcorner = geo_reference.get_yllcorner()
-        zone = geo_reference.get_zone()
-    else:
-        zone = origin[0]
-        xllcorner = origin[1]
-        yllcorner = origin[2]
-
-    # FIXME: Refactor using code from Interpolation_function.statistics
-    # (in interpolate.py)
-    # Something like print swwstats(swwname)
-    if verbose:
-        log.critical('------------------------------------------------')
-        log.critical('Statistics of SWW file:')
-        log.critical('  Name: %s' % swwfile)
-        log.critical('  Reference:')
-        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
-        if type(reduction) is not types.BuiltinFunctionType:
-            log.critical('    Time: %f' % times)
-        else:
-            log.critical('    Start time: %f' % fid.starttime[0])
-        log.critical('  Extent:')
-        log.critical('    x [m] in [%f, %f], len(x) == %d'
-                     %(num.min(x), num.max(x), len(x.flat)))
-        log.critical('    y [m] in [%f, %f], len(y) == %d'
-                     % (num.min(y), num.max(y), len(y.flat)))
-        if type(reduction) is not types.BuiltinFunctionType:
-            log.critical('    t [s] = %f, len(t) == %d' % (times, 1))
-        else:
-            log.critical('    t [s] in [%f, %f], len(t) == %d'
-                         % (min(times), max(times), len(times)))
-        log.critical('  Quantities [SI units]:')
-        
-        # Comment out for reduced memory consumption
-        for name in ['stage', 'xmomentum', 'ymomentum']:
-            q = fid.variables[name][:].flatten()
-            if type(reduction) is not types.BuiltinFunctionType:
-                q = q[reduction*len(x):(reduction+1)*len(x)]
-            if verbose: log.critical('    %s in [%f, %f]'
-                                     % (name, min(q), max(q)))
-        for name in ['elevation']:
-            q = fid.variables[name][:].flatten()
-            if verbose: log.critical('    %s in [%f, %f]'
-                                     % (name, min(q), max(q)))
-
-    # Get the variables in the supplied expression.
-    # This may throw a SyntaxError exception.
-    var_list = get_vars_in_expression(quantity)
-
-    # Check that we have the required variables in the SWW file.
-    missing_vars = []
-    for name in var_list:
-        try:
-            _ = fid.variables[name]
-        except:
-            missing_vars.append(name)
-    if missing_vars:
-        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
-               % (quantity, swwfile))
-        raise Exception, msg
-
-    # Create result array and start filling, block by block.
-    result = num.zeros(number_of_points, num.float)
-
-    for start_slice in xrange(0, number_of_points, block_size):
-        # Limit slice size to array end if at last block
-        end_slice = min(start_slice + block_size, number_of_points)
-        
-        # Get slices of all required variables
-        q_dict = {}
-        for name in var_list:
-            # check if variable has time axis
-            if len(fid.variables[name].shape) == 2:
-                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
-            else:       # no time axis
-                q_dict[name] = fid.variables[name][start_slice:end_slice]
-
-        # Evaluate expression with quantities found in SWW file
-        res = apply_expression_to_dictionary(quantity, q_dict)
-
-        if len(res.shape) == 2:
-            new_res = num.zeros(res.shape[1], num.float)
-            for k in xrange(res.shape[1]):
-                if type(reduction) is not types.BuiltinFunctionType:
-                    new_res[k] = res[reduction,k]
-                else:
-                    new_res[k] = reduction(res[:,k])
-            res = new_res
-
-        result[start_slice:end_slice] = res
-                                    
-    # Post condition: Now q has dimension: number_of_points
-    assert len(result.shape) == 1
-    assert result.shape[0] == number_of_points
-
-    if verbose:
-        log.critical('Processed values for %s are in [%f, %f]'
-                     % (quantity, min(result), max(result)))
-
-    # Create grid and update xll/yll corner and x,y
-    # Relative extent
-    if easting_min is None:
-        xmin = min(x)
-    else:
-        xmin = easting_min - xllcorner
-
-    if easting_max is None:
-        xmax = max(x)
-    else:
-        xmax = easting_max - xllcorner
-
-    if northing_min is None:
-        ymin = min(y)
-    else:
-        ymin = northing_min - yllcorner
-
-    if northing_max is None:
-        ymax = max(y)
-    else:
-        ymax = northing_max - yllcorner
-
-    msg = 'xmax must be greater than or equal to xmin.\n'
-    msg += 'I got xmin = %f, xmax = %f' %(xmin, xmax)
-    assert xmax >= xmin, msg
-
-    msg = 'ymax must be greater than or equal to xmin.\n'
-    msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
-    assert ymax >= ymin, msg
-
-    if verbose: log.critical('Creating grid')
-    ncols = int((xmax-xmin)/cellsize) + 1
-    nrows = int((ymax-ymin)/cellsize) + 1
-
-    # New absolute reference and coordinates
-    newxllcorner = xmin + xllcorner
-    newyllcorner = ymin + yllcorner
-
-    x = x + xllcorner - newxllcorner
-    y = y + yllcorner - newyllcorner
-
-    vertex_points = num.concatenate ((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
-    assert len(vertex_points.shape) == 2
-
-    grid_points = num.zeros ((ncols*nrows, 2), num.float)
-
-    for i in xrange(nrows):
-        if format.lower() == 'asc':
-            yg = i * cellsize
-        else:
-            # this will flip the order of the y values for ers
-            yg = (nrows-i) * cellsize
-
-        for j in xrange(ncols):
-            xg = j * cellsize
-            k = i*ncols + j
-
-            grid_points[k, 0] = xg
-            grid_points[k, 1] = yg
-
-    # Interpolate
-    from anuga.fit_interpolate.interpolate import Interpolate
-
-    # Remove loners from vertex_points, volumes here
-    vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
-    # export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
-    interp = Interpolate(vertex_points, volumes, verbose = verbose)
-
-    # Interpolate using quantity values
-    if verbose: log.critical('Interpolating')
-    grid_values = interp.interpolate(result, grid_points).flatten()
-
-    if verbose:
-        log.critical('Interpolated values are in [%f, %f]'
-                     % (num.min(grid_values), num.max(grid_values)))
-
-    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
-    P = interp.mesh.get_boundary_polygon()
-    outside_indices = outside_polygon(grid_points, P, closed=True)
-
-    for i in outside_indices:
-        grid_values[i] = NODATA_value
-
-    if format.lower() == 'ers':
-        # setup ERS header information
-        grid_values = num.reshape(grid_values, (nrows, ncols))
-        header = {}
-        header['datum'] = '"' + datum + '"'
-        # FIXME The use of hardwired UTM and zone number needs to be made optional
-        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
-        header['projection'] = '"UTM-' + str(zone) + '"'
-        header['coordinatetype'] = 'EN'
-        if header['coordinatetype'] == 'LL':
-            header['longitude'] = str(newxllcorner)
-            header['latitude'] = str(newyllcorner)
-        elif header['coordinatetype'] == 'EN':
-            header['eastings'] = str(newxllcorner)
-            header['northings'] = str(newyllcorner)
-        header['nullcellvalue'] = str(NODATA_value)
-        header['xdimension'] = str(cellsize)
-        header['ydimension'] = str(cellsize)
-        header['value'] = '"' + quantity + '"'
-        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
-
-        #Write
-        if verbose: log.critical('Writing %s' % demfile)
-
-        import ermapper_grids
-
-        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
-
-        fid.close()
-    else:
-        #Write to Ascii format
-        #Write prj file
-        prjfile = basename_out + '.prj'
-
-        if verbose: log.critical('Writing %s' % prjfile)
-        prjid = open(prjfile, 'w')
-        prjid.write('Projection    %s\n' %'UTM')
-        prjid.write('Zone          %d\n' %zone)
-        prjid.write('Datum         %s\n' %datum)
-        prjid.write('Zunits        NO\n')
-        prjid.write('Units         METERS\n')
-        prjid.write('Spheroid      %s\n' %datum)
-        prjid.write('Xshift        %d\n' %false_easting)
-        prjid.write('Yshift        %d\n' %false_northing)
-        prjid.write('Parameters\n')
-        prjid.close()
-
-        if verbose: log.critical('Writing %s' % demfile)
-
-        ascid = open(demfile, 'w')
-
-        ascid.write('ncols         %d\n' %ncols)
-        ascid.write('nrows         %d\n' %nrows)
-        ascid.write('xllcorner     %d\n' %newxllcorner)
-        ascid.write('yllcorner     %d\n' %newyllcorner)
-        ascid.write('cellsize      %f\n' %cellsize)
-        ascid.write('NODATA_value  %d\n' %NODATA_value)
-
-        #Get bounding polygon from mesh
-        #P = interp.mesh.get_boundary_polygon()
-        #inside_indices = inside_polygon(grid_points, P)
-
-        # change printoptions so that a long string of zeros in not
-        # summarized as [0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0]
-        #printoptions = num.get_printoptions()
-        #num.set_printoptions(threshold=sys.maxint)
-
-        format = '%.'+'%g' % number_of_decimal_places +'e'
-        for i in range(nrows):
-            if verbose and i % ((nrows+10)/10) == 0:
-                log.critical('Doing row %d of %d' % (i, nrows))
-
-            base_index = (nrows-i-1)*ncols
-
-            slice = grid_values[base_index:base_index+ncols]
-
-            num.savetxt(ascid, slice.reshape(1,ncols), format, ' ' )
-            
-        
-        #Close
-        ascid.close()
-        fid.close()
-
-        return basename_out
-
-################################################################################
-# Obsolete functions - maintained for backwards compatibility
-################################################################################
-
-##
-# @brief 
-# @param basename_in 
-# @param basename_out 
-# @param quantity 
-# @param timestep 
-# @param reduction 
-# @param cellsize 
-# @param verbose 
-# @param origin 
-# @note OBSOLETE - use sww2dem() instead.
-def sww2asc(basename_in, basename_out = None,
-            quantity = None,
-            reduction = None,
-            cellsize = 10,
-            verbose = False,
-            origin = None):
-    
-    log.critical('sww2asc will soon be obsoleted - please use sww2dem')
-    sww2dem(basename_in,
-            basename_out = basename_out,
-            quantity = quantity,
-            reduction = reduction,
-            cellsize = cellsize,
-            number_of_decimal_places = number_of_decimal_places,
-            verbose = verbose,
-            origin = origin,
-            datum = 'WGS84',
-            format = 'asc')
-
-
-##
-# @brief 
-# @param basename_in
-# @param basename_out
-# @param quantity
-# @param timestep
-# @param reduction
-# @param cellsize
-# @param verbose
-# @param origin
-# @param datum
-# @note OBSOLETE - use sww2dem() instead.
-def sww2ers(basename_in, basename_out=None,
-            quantity=None,
-            reduction=None,
-            cellsize=10,
-            verbose=False,
-            origin=None,
-            datum='WGS84'):
-    log.critical('sww2ers will soon be obsoleted - please use sww2dem')
-    sww2dem(basename_in,
-            basename_out=basename_out,
-            quantity=quantity,
-            reduction=reduction,
-            cellsize=cellsize,
-            number_of_decimal_places=number_of_decimal_places,
-            verbose=verbose,
-            origin=origin,
-            datum=datum,
-            format='ers')
 
 ################################################################################
@@ -3544,124 +3026,4 @@
 ################################################################################
 
-##
-# @brief Store screen output and errors to a file.
-# @param dir_name Path to directory for output files (default '.').
-# @param myid 
-# @param numprocs 
-# @param extra_info 
-# @param verbose True if this function is to be verbose.
-def start_screen_catcher(dir_name=None, myid='', numprocs='', extra_info='',
-                         verbose=True):
-    """
-    Used to store screen output and errors to file, if run on multiple
-    processes each processor will have its own output and error file.
-
-    extra_info - is used as a string that can identify outputs with another
-    string eg. '_other'
-
-    FIXME: Would be good if you could suppress all the screen output and
-    only save it to file... however it seems a bit tricky as this capture
-    techique response to sys.stdout and by this time it is already printed out.
-    """
-
-    import sys
-
-    if dir_name == None:
-        dir_name = getcwd()
-
-    if access(dir_name, W_OK) == 0:
-        if verbose: log.critical('Making directory %s' % dir_name)
-        mkdir (dir_name,0777)
-
-    if myid != '':
-        myid = '_' + str(myid)
-    if numprocs != '':
-        numprocs = '_' + str(numprocs)
-    if extra_info != '':
-        extra_info = '_' + str(extra_info)
-
-    screen_output_name = join(dir_name, "screen_output%s%s%s.txt" %
-                                        (myid, numprocs, extra_info))
-    screen_error_name = join(dir_name, "screen_error%s%s%s.txt" %
-                                       (myid, numprocs, extra_info))
-
-    if verbose: log.critical('Starting ScreenCatcher, all output will be '
-                             'stored in %s' % screen_output_name)
-
-    # used to catch screen output to file
-    sys.stdout = Screen_Catcher(screen_output_name)
-    sys.stderr = Screen_Catcher(screen_error_name)
-
-
-##
-# @brief A class to catch stdout and stderr and write to files.
-class Screen_Catcher:
-    """this simply catches the screen output and stores it to file defined by
-    start_screen_catcher (above)
-    """
-
-    ##
-    # @brief Initialize this instance of Screen_Catcher.
-    # @param filename The path to the file to write to.
-    def __init__(self, filename):
-        self.filename = filename
-        if exists(self.filename)is True:
-            log.critical('Old existing file "%s" has been deleted'
-                         % self.filename)
-            remove(self.filename)
-
-    ##
-    # @brief Write output to the file.
-    # @param stuff The string to write.
-    def write(self, stuff):
-        fid = open(self.filename, 'a')
-        fid.write(stuff)
-        fid.close()
-
-
-##
-# @brief Copy a file to a directory, and optionally append another file to it.
-# @param dir_name Target directory.
-# @param filename Path to file to copy to directory 'dir_name'.
-# @param filename2 Optional path to file to append to copied file.
-# @param verbose True if this function is to be verbose.
-# @note Allow filenames to be either a string or sequence of strings.
-def copy_code_files(dir_name, filename1, filename2=None, verbose=False):
-    """Copies "filename1" and "filename2" to "dir_name".
-
-    Each 'filename' may be a string or list of filename strings.
-
-    Filenames must be absolute pathnames
-    """
-
-    ##
-    # @brief copies a file or sequence to destination directory.
-    # @param dest The destination directory to copy to.
-    # @param file A filename string or sequence of filename strings.
-    def copy_file_or_sequence(dest, file):
-        if hasattr(file, '__iter__'):
-            for f in file:
-                shutil.copy(f, dir_name)
-                if verbose:
-                    log.critical('File %s copied' % f)
-        else:
-            shutil.copy(file, dir_name)
-            if verbose:
-                log.critical('File %s copied' % file)
-
-    # check we have a destination directory, create if necessary
-    if not os.path.isdir(dir_name):
-        if verbose:
-            log.critical('Make directory %s' % dir_name)
-        mkdir(dir_name, 0777)
-
-    if verbose:
-        log.critical('Output directory: %s' % dir_name)        
-
-    copy_file_or_sequence(dir_name, filename1)
-
-    if not filename2 is None:
-        copy_file_or_sequence(dir_name, filename2)
-
 
 ##
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7743)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7744)
@@ -20,9 +20,10 @@
 from anuga.file_conversion.file_conversion import tsh2sww, \
                         asc_csiro2sww, pmesh_to_domain_instance
-                        
+from anuga.coordinate_transforms.geo_reference import Geo_reference                        
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
 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.utilities.file_utils import del_dir, load_csv_as_dict, \
+                                        load_csv_as_array
 from anuga.anuga_exceptions import ANUGAError
 from anuga.utilities.numerical_tools import ensure_numeric, mean
@@ -866,202 +867,4 @@
 
 
-    def test_sww2dem_asc_elevation_depth(self):
-        """test_sww2dem_asc_elevation_depth
-        
-        Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        
-        Also check geo_reference is correct
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        self.domain.set_name('datatest')
-
-        prjfile = self.domain.get_name() + '_elevation.prj'
-        ascfile = self.domain.get_name() + '_elevation.asc'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 1.0)
-
-        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-        
-        # Check georeferencig: zone, xllcorner and yllcorner
-        assert fid.zone[0] == 56
-        assert fid.xllcorner[0] == 308500
-        assert fid.yllcorner[0] == 6189000
-                
-
-        fid.close()
-
-        #Export to ascii/prj files
-        sww2dem(self.domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                verbose = self.verbose,
-                format = 'asc')
-
-        #Check prj (meta data)
-        prjid = open(prjfile)
-        lines = prjid.readlines()
-        prjid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'projection'
-        assert L[1].strip().lower() == 'utm'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'zone'
-        assert L[1].strip().lower() == '56'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'datum'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'zunits'
-        assert L[1].strip().lower() == 'no'
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'units'
-        assert L[1].strip().lower() == 'meters'
-
-        L = lines[5].strip().split()
-        assert L[0].strip().lower() == 'spheroid'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[6].strip().split()
-        assert L[0].strip().lower() == 'xshift'
-        assert L[1].strip().lower() == '500000'
-
-        L = lines[7].strip().split()
-        assert L[0].strip().lower() == 'yshift'
-        assert L[1].strip().lower() == '10000000'
-
-        L = lines[8].strip().split()
-        assert L[0].strip().lower() == 'parameters'
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), -i*cellsize - y)
-                
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-
-        #Export to ascii/prj files
-        sww2dem(self.domain.get_name(),
-                quantity = 'depth',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                verbose = self.verbose,
-                format = 'asc')
-        
-        #Check asc file
-        ascfile = self.domain.get_name() + '_depth.asc'
-        prjfile = self.domain.get_name() + '_depth.prj'
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
-
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
     def test_export_grid(self):
         """
@@ -1556,1125 +1359,4 @@
         os.remove(ascfile)
         os.remove(swwfile2)
-
-    def test_sww2dem_larger(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView. Here:
-
-        ncols         11
-        nrows         11
-        xllcorner     308500
-        yllcorner     6189000
-        cellsize      10.000000
-        NODATA_value  -9999
-        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
-         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
-         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
-         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
-         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
-         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
-         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
-         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
-         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
-         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
-           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
-
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-
-        from mesh_factory import rectangular
-
-        #Create basic mesh (100m x 100m)
-        points, vertices, boundary = rectangular(2, 2, 100, 100)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-
-        domain.set_name('datatest')
-
-        prjfile = domain.get_name() + '_elevation.prj'
-        ascfile = domain.get_name() + '_elevation.asc'
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-        domain.geo_reference = Geo_reference(56, 308500, 6189000)
-
-        #
-        domain.set_quantity('elevation', lambda x,y: -x-y)
-        domain.set_quantity('stage', 0)
-
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
-
-
-        #
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        
-        domain.tight_slope_limiters = 1
-        domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 10  #10m grid
-
-
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ascii/prj files
-        sww2dem(domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                verbose = self.verbose,
-                format = 'asc',
-                block_size=2)
-
-
-        #Check prj (meta data)
-        prjid = open(prjfile)
-        lines = prjid.readlines()
-        prjid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'projection'
-        assert L[1].strip().lower() == 'utm'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'zone'
-        assert L[1].strip().lower() == '56'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'datum'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'zunits'
-        assert L[1].strip().lower() == 'no'
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'units'
-        assert L[1].strip().lower() == 'meters'
-
-        L = lines[5].strip().split()
-        assert L[0].strip().lower() == 'spheroid'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[6].strip().split()
-        assert L[0].strip().lower() == 'xshift'
-        assert L[1].strip().lower() == '500000'
-
-        L = lines[7].strip().split()
-        assert L[0].strip().lower() == 'yshift'
-        assert L[1].strip().lower() == '10000000'
-
-        L = lines[8].strip().split()
-        assert L[0].strip().lower() == 'parameters'
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '11'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '11'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
-        for i, line in enumerate(lines[6:]):
-            for j, value in enumerate( line.split() ):
-                assert num.allclose(float(value), -(10-i+j)*cellsize,
-                                    atol=1.0e-12, rtol=1.0e-12)
-
-                # Note: Equality can be obtained in this case,
-                # but it is better to use allclose.
-                #assert float(value) == -(10-i+j)*cellsize
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-
-    def test_sww2dem_larger_zero(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. Arcview. This example has rows with a
-        large number of zeros
-
-        ncols         2001
-        nrows         2
-        xllcorner     308500
-        yllcorner     6189000
-        cellsize      1.000000
-        NODATA_value  -9999
-        0.0 ....
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-
-        from mesh_factory import rectangular_cross
-
-        #Create basic mesh (100m x 100m)
-        points, vertices, boundary = rectangular_cross(2000, 1, 2000.0, 1.0)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 1
-
-        domain.set_name('datatest')
-
-        prjfile = domain.get_name() + '_elevation.prj'
-        ascfile = domain.get_name() + '_elevation.asc'
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-        domain.geo_reference = Geo_reference(56, 308500, 6189000)
-
-        #
-        domain.set_quantity('elevation', 0)
-        domain.set_quantity('stage', 0)
-
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
-
-
-        #
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        
-        domain.tight_slope_limiters = 1
-        domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 1.0  #0.1 grid
-
-
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ascii/prj files
-        sww2dem(domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                verbose = self.verbose,
-                format = 'asc',
-                block_size=2)
-
-
-        #Check prj (meta data)
-        prjid = open(prjfile)
-        lines = prjid.readlines()
-        prjid.close()
-
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'projection'
-        assert L[1].strip().lower() == 'utm'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'zone'
-        assert L[1].strip().lower() == '56'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'datum'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'zunits'
-        assert L[1].strip().lower() == 'no'
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'units'
-        assert L[1].strip().lower() == 'meters'
-
-        L = lines[5].strip().split()
-        assert L[0].strip().lower() == 'spheroid'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[6].strip().split()
-        assert L[0].strip().lower() == 'xshift'
-        assert L[1].strip().lower() == '500000'
-
-        L = lines[7].strip().split()
-        assert L[0].strip().lower() == 'yshift'
-        assert L[1].strip().lower() == '10000000'
-
-        L = lines[8].strip().split()
-        assert L[0].strip().lower() == 'parameters'
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '2001'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '2'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
-        for i, line in enumerate(lines[6:]):
-            for j, value in enumerate( line.split() ):
-                assert num.allclose(float(value), 0.0,
-                                    atol=1.0e-12, rtol=1.0e-12)
-
-                # Note: Equality can be obtained in this case,
-                # but it is better to use allclose.
-                #assert float(value) == -(10-i+j)*cellsize
-
-
-        fid.close()
-
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-
-
-    def test_sww2dem_boundingbox(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView.
-        This will test that mesh can be restricted by bounding box
-
-        Original extent is 100m x 100m:
-
-        Eastings:   308500 -  308600
-        Northings: 6189000 - 6189100
-
-        Bounding box changes this to the 50m x 50m square defined by
-
-        Eastings:   308530 -  308570
-        Northings: 6189050 - 6189100
-
-        The cropped values should be
-
-         -130 -140 -150 -160 -170
-         -120 -130 -140 -150 -160
-         -110 -120 -130 -140 -150
-         -100 -110 -120 -130 -140
-          -90 -100 -110 -120 -130
-          -80  -90 -100 -110 -120
-
-        and the new lower reference point should be
-        Eastings:   308530
-        Northings: 6189050
-
-        Original dataset is the same as in test_sww2dem_larger()
-
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-
-        from mesh_factory import rectangular
-
-        #Create basic mesh (100m x 100m)
-        points, vertices, boundary = rectangular(2, 2, 100, 100)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-
-        domain.set_name('datatest')
-
-        prjfile = domain.get_name() + '_elevation.prj'
-        ascfile = domain.get_name() + '_elevation.asc'
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-        domain.geo_reference = Geo_reference(56, 308500, 6189000)
-
-        #
-        domain.set_quantity('elevation', lambda x,y: -x-y)
-        domain.set_quantity('stage', 0)
-
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
-
-
-        #
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #domain.tight_slope_limiters = 1
-        domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 10  #10m grid
-
-
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        # Export to ascii/prj files
-        sww2dem(domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                easting_min = 308530,
-                easting_max = 308570,
-                northing_min = 6189050,
-                northing_max = 6189100,
-                verbose = self.verbose,
-                format = 'asc')
-
-        fid.close()
-
-
-        # Check prj (meta data)
-        prjid = open(prjfile)
-        lines = prjid.readlines()
-        prjid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'projection'
-        assert L[1].strip().lower() == 'utm'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'zone'
-        assert L[1].strip().lower() == '56'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'datum'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'zunits'
-        assert L[1].strip().lower() == 'no'
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'units'
-        assert L[1].strip().lower() == 'meters'
-
-        L = lines[5].strip().split()
-        assert L[0].strip().lower() == 'spheroid'
-        assert L[1].strip().lower() == 'wgs84'
-
-        L = lines[6].strip().split()
-        assert L[0].strip().lower() == 'xshift'
-        assert L[1].strip().lower() == '500000'
-
-        L = lines[7].strip().split()
-        assert L[0].strip().lower() == 'yshift'
-        assert L[1].strip().lower() == '10000000'
-
-        L = lines[8].strip().split()
-        assert L[0].strip().lower() == 'parameters'
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '6'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308530)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189050)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values
-        for i, line in enumerate(lines[6:]):
-            for j, value in enumerate( line.split() ):
-                #assert float(value) == -(10-i+j)*cellsize
-                assert float(value) == -(10-i+j+3)*cellsize
-
-
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-
-    def test_sww2dem_asc_stage_reduction(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-
-        This tests the reduction of quantity stage using min
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-        self.domain.set_name('datatest')
-
-        prjfile = self.domain.get_name() + '_stage.prj'
-        ascfile = self.domain.get_name() + '_stage.asc'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ascii/prj files
-        sww2dem(self.domain.get_name(),
-                quantity = 'stage',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                reduction = min,
-                format = 'asc',
-                verbose=self.verbose)
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-
-        #Check grid values (where applicable)
-        for j in range(5):
-            if j%2 == 0:
-                L = lines[6+j].strip().split()
-                jj = 4-j
-                for i in range(5):
-                    if i%2 == 0:
-                        index = jj/2 + i/2*3
-                        val0 = stage[0,index]
-                        val1 = stage[1,index]
-
-                        #print i, j, index, ':', L[i], val0, val1
-                        assert num.allclose(float(L[i]), min(val0, val1))
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-    def test_sww2dem_asc_stage_time(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-
-        This tests the reduction of quantity stage using min
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-        self.domain.set_name('datatest')
-
-        prjfile = self.domain.get_name() + '_stage.prj'
-        ascfile = self.domain.get_name() + '_stage.asc'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-        #Export to ascii/prj files
-        sww2dem(self.domain.get_name(),
-                quantity = 'stage',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                reduction = 1,
-                format = 'asc',
-                verbose=self.verbose)
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values (where applicable)
-        for j in range(5):
-            if j%2 == 0:
-                L = lines[6+j].strip().split()
-                jj = 4-j
-                for i in range(5):
-                    if i%2 == 0:
-                        index = jj/2 + i/2*3
-                        
-                        val = stage[1,index]
-                   
-                        assert num.allclose(float(L[i]), val)
-
-        fid.close()
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-    def test_sww2dem_asc_derived_quantity(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-
-        This tests the use of derived quantities
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-        self.domain.set_name('datatest')
-
-        prjfile = self.domain.get_name() + '_depth.prj'
-        ascfile = self.domain.get_name() + '_depth.asc'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 0.0)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ascii/prj files
-        sww2dem(self.domain.get_name(),
-                basename_out = 'datatest_depth',
-                quantity = 'stage - elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                reduction = min,
-                format = 'asc',
-                verbose = self.verbose)
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-
-        #Check grid values (where applicable)
-        for j in range(5):
-            if j%2 == 0:
-                L = lines[6+j].strip().split()
-                jj = 4-j
-                for i in range(5):
-                    if i%2 == 0:
-                        index = jj/2 + i/2*3
-                        val0 = stage[0,index] - z[index]
-                        val1 = stage[1,index] - z[index]
-
-                        #print i, j, index, ':', L[i], val0, val1
-                        assert num.allclose(float(L[i]), min(val0, val1))
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-
-
-
-    def test_sww2dem_asc_missing_points(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-
-        This test includes the writing of missing values
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup mesh not coinciding with rectangle.
-        #This will cause missing values to occur in gridded data
-
-
-        points = [                        [1.0, 1.0],
-                              [0.5, 0.5], [1.0, 0.5],
-                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
-
-        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
-
-        #Create shallow water domain
-        domain = Domain(points, vertices)
-        domain.default_order=2
-
-
-        #Set some field values
-        domain.set_quantity('elevation', lambda x,y: -x-y)
-        domain.set_quantity('friction', 0.03)
-
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'exterior': B} )
-
-
-        ######################
-        #Initial condition - with jumps
-
-        bed = domain.quantities['elevation'].vertex_values
-        stage = num.zeros(bed.shape, num.float)
-
-        h = 0.3
-        for i in range(stage.shape[0]):
-            if i % 2 == 0:
-                stage[i,:] = bed[i,:] + h
-            else:
-                stage[i,:] = bed[i,:]
-
-        domain.set_quantity('stage', stage)
-        domain.distribute_to_vertices_and_edges()
-
-        domain.set_name('datatest')
-
-        prjfile = domain.get_name() + '_elevation.prj'
-        ascfile = domain.get_name() + '_elevation.asc'
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(swwfile, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
-
-        #Export to ascii/prj files
-        sww2dem(domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                verbose = self.verbose,
-                format = 'asc')
-
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[0].strip().split()
-        assert L[0].strip().lower() == 'ncols'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[1].strip().split()
-        assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        L = lines[4].strip().split()
-        assert L[0].strip().lower() == 'cellsize'
-        assert num.allclose(float(L[1].strip().lower()), cellsize)
-
-        L = lines[5].strip().split()
-        assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            assert len(L) == 5
-            y = (4-j) * cellsize
-
-            for i in range(5):
-                #print i
-                if i+j >= 4:
-                    assert num.allclose(float(L[i]), -i*cellsize - y)
-                else:
-                    #Missing values
-                    assert num.allclose(float(L[i]), -9999)
-
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-    def test_sww2ers_simple(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-
-        NODATA_value = 1758323
-
-        #Setup
-        self.domain.set_name('datatest')
-
-        headerfile = self.domain.get_name() + '.ers'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ers files
-        sww2dem(self.domain.get_name(),
-                quantity = 'elevation',
-                cellsize = cellsize,
-                number_of_decimal_places = 9,
-                NODATA_value = NODATA_value,
-                verbose = self.verbose,
-                format = 'ers')
-
-        #Check header data
-        from ermapper_grids import read_ermapper_header, read_ermapper_data
-
-        header = read_ermapper_header(self.domain.get_name() + '_elevation.ers')
-        #print header
-        assert header['projection'].lower() == '"utm-56"'
-        assert header['datum'].lower() == '"wgs84"'
-        assert header['units'].lower() == '"meters"'
-        assert header['value'].lower() == '"elevation"'
-        assert header['xdimension'] == '0.25'
-        assert header['ydimension'] == '0.25'
-        assert float(header['eastings']) == 308500.0   #xllcorner
-        assert float(header['northings']) == 6189000.0 #yllcorner
-        assert int(header['nroflines']) == 5
-        assert int(header['nrofcellsperline']) == 5
-        assert int(header['nullcellvalue']) == NODATA_value
-        #FIXME - there is more in the header
-
-
-        #Check grid data
-        grid = read_ermapper_data(self.domain.get_name() + '_elevation')
-
-        #FIXME (Ole): Why is this the desired reference grid for -x-y?
-        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
-                    -1,    -1.25, -1.5,  -1.75, -2.0,
-                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
-                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
-                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
-
-
-        #print grid
-        assert num.allclose(grid, ref_grid)
-
-        fid.close()
-
-        #Cleanup
-        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
-        #Done (Ole) - it was because sww2ers didn't close it's sww file
-        os.remove(sww.filename)
-        os.remove(self.domain.get_name() + '_elevation')
-        os.remove(self.domain.get_name() + '_elevation.ers')
 
 
@@ -3347,127 +2029,5 @@
         os.remove(root + '_100.dem')
 
-    def xxxtestz_sww2ers_real(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-
-        # the memory optimised least squares
-        #  cellsize = 20,   # this one seems to hang
-        #  cellsize = 200000, # Ran 1 test in 269.703s
-                                #Ran 1 test in 267.344s
-        #  cellsize = 20000,  # Ran 1 test in 460.922s
-        #  cellsize = 2000   #Ran 1 test in 5340.250s
-        #  cellsize = 200   #this one seems to hang, building matirx A
-
-        # not optimised
-        # seems to hang
-        #  cellsize = 2000   # Ran 1 test in 5334.563s
-        #Export to ascii/prj files
-        sww2dem('karratha_100m',
-                quantity = 'depth',
-                cellsize = 200000,
-                number_of_decimal_places = 9,
-                verbose = True)
-
-    def test_read_asc(self):
-        """Test conversion from dem in ascii format to native NetCDF format
-        """
-
-        import time, os
-
-        from file_conversion import _read_asc
-        #Write test asc file
-        filename = tempfile.mktemp(".000")
-        fid = open(filename, 'w')
-        fid.write("""ncols         7
-nrows         4
-xllcorner     2000.5
-yllcorner     3000.5
-cellsize      25
-NODATA_value  -9999
-    97.921    99.285   125.588   180.830   258.645   342.872   415.836
-   473.157   514.391   553.893   607.120   678.125   777.283   883.038
-   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
-   502.645   516.230   504.739   450.604   388.500   338.097   514.980
-""")
-        fid.close()
-        bath_metadata, grid = _read_asc(filename, verbose=self.verbose)
-        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
-        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
-        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
-        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
-        self.failUnless(grid[0][0]  == 97.921,  'Failed')
-        self.failUnless(grid[3][6]  == 514.980,  'Failed')
-
-        os.remove(filename)
-
-
-    #### TESTS FOR URS 2 SWW  ###     
-    
-    def create_mux(self, points_num=None):
-        # write all the mux stuff.
-        time_step_count = 3
-        time_step = 0.5
-        
-        longitudes = [150.66667, 150.83334, 151., 151.16667]
-        latitudes = [-34.5, -34.33333, -34.16667, -34]
-
-        if points_num == None:
-            points_num = len(longitudes) * len(latitudes)
-
-        lonlatdeps = []
-        quantities = ['HA','UA','VA']
-        mux_names = [WAVEHEIGHT_MUX_LABEL,
-                     EAST_VELOCITY_LABEL,
-                     NORTH_VELOCITY_LABEL]
-        quantities_init = [[],[],[]]
-        # urs binary is latitude fastest
-        for i,lon in enumerate(longitudes):
-            for j,lat in enumerate(latitudes):
-                _ , e, n = redfearn(lat, lon)
-                lonlatdeps.append([lon, lat, n])
-                quantities_init[0].append(e) # HA
-                quantities_init[1].append(n ) # UA
-                quantities_init[2].append(e) # VA
-        #print "lonlatdeps",lonlatdeps
-
-        file_handle, base_name = tempfile.mkstemp("")
-        os.close(file_handle)
-        os.remove(base_name)
-        
-        files = []        
-        for i,q in enumerate(quantities): 
-            quantities_init[i] = ensure_numeric(quantities_init[i])
-            #print "HA_init", HA_init
-            q_time = num.zeros((time_step_count, points_num), num.float)
-            for time in range(time_step_count):
-                q_time[time,:] = quantities_init[i] #* time * 4
-            
-            #Write C files
-            columns = 3 # long, lat , depth
-            file = base_name + mux_names[i]
-            #print "base_name file",file 
-            f = open(file, 'wb')
-            files.append(file)
-            f.write(pack('i',points_num))
-            f.write(pack('i',time_step_count))
-            f.write(pack('f',time_step))
-
-            #write lat/long info
-            for lonlatdep in lonlatdeps:
-                for float in lonlatdep:
-                    f.write(pack('f',float))
-                    
-            # Write quantity info
-            for time in  range(time_step_count):
-                for point_i in range(points_num):
-                    f.write(pack('f',q_time[time,point_i]))
-                    #print " mux_names[i]", mux_names[i] 
-                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
-            f.close()
-        return base_name, files
-    
+
     def write_mux(self, lat_long_points, time_step_count, time_step,
                   depth=None, ha=None, ua=None, va=None):
Index: /anuga_core/source/anuga/shallow_water/test_forcing.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_forcing.py	(revision 7743)
+++ /anuga_core/source/anuga/shallow_water/test_forcing.py	(revision 7744)
@@ -6,4 +6,5 @@
 from boundaries import Reflective_boundary
 from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.file_conversion.file_conversion import timefile2netcdf
 from forcing import *
 
@@ -246,7 +247,4 @@
 
         fid.close()
-
-        # Convert ASCII file to NetCDF (Which is what we really like!)
-        from file_conversion import timefile2netcdf
 
         timefile2netcdf(filename)
@@ -338,7 +336,4 @@
 
         fid.close()
-
-        # Convert ASCII file to NetCDF (Which is what we really like!)
-        from file_conversion import timefile2netcdf
 
         timefile2netcdf(filename, time_as_seconds=True)
Index: /anuga_core/source/anuga/utilities/file_utils.py
===================================================================
--- /anuga_core/source/anuga/utilities/file_utils.py	(revision 7743)
+++ /anuga_core/source/anuga/utilities/file_utils.py	(revision 7744)
@@ -6,4 +6,5 @@
 import os, sys
 import csv
+import numpy as num
 
 ##
@@ -362,2 +363,30 @@
 
     return attribute_dic, title_index_dic
+
+
+          
+##
+# @brief Convert CSV file to a dictionary of arrays.
+# @param file_name The path to the file to read.
+def load_csv_as_array(file_name):
+    """
+    Convert CSV files of the form:
+
+    time, discharge, velocity
+    0.0,  1.2,       0.0
+    0.1,  3.2,       1.1
+    ...
+
+    to a dictionary of numeric arrays.
+
+
+    See underlying function load_csv_as_dict for more details.
+    """
+
+    X, _ = load_csv_as_dict(file_name)
+
+    Y = {}
+    for key in X.keys():
+        Y[key] = num.array([float(x) for x in X[key]])
+
+    return Y
