Index: /trunk/anuga_core/documentation/user_manual/demos/channel1.py
===================================================================
--- /trunk/anuga_core/documentation/user_manual/demos/channel1.py	(revision 7769)
+++ /trunk/anuga_core/documentation/user_manual/demos/channel1.py	(revision 7770)
@@ -7,8 +7,11 @@
 # Import necessary modules
 #------------------------------------------------------------------------------
+
+# Import standard shallow water domain and standard boundaries.
+import anuga
+
+# Import rectangular cross mesh generator
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Dirichlet_boundary
+
 
 #------------------------------------------------------------------------------
@@ -18,5 +21,5 @@
                                                len1=10.0, len2=5.0) # Mesh
 
-domain = Domain(points, vertices, boundary)  # Create domain
+domain = anuga.Domain(points, vertices, boundary)  # Create domain
 domain.set_name('channel1')                  # Output name
 
@@ -35,6 +38,6 @@
 # Setup boundary conditions
 #------------------------------------------------------------------------------
-Bi = Dirichlet_boundary([0.4, 0, 0])         # Inflow
-Br = Reflective_boundary(domain)             # Solid reflective wall
+Bi = anuga.Dirichlet_boundary([0.4, 0, 0])         # Inflow
+Br = anuga.Reflective_boundary(domain)             # Solid reflective wall
 
 domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
Index: /trunk/anuga_core/source/anuga/file/sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/sww.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file/sww.py	(revision 7770)
@@ -1196,4 +1196,92 @@
 
 ##
+# @brief Get mesh and quantity data from an SWW file.
+# @param filename Path to data file to read.
+# @param quantities UNUSED!
+# @param verbose True if this function is to be verbose.
+# @return (mesh, quantities, time) where mesh is the mesh data, quantities is
+#         a dictionary of {name: value}, and time is the time vector.
+# @note Quantities extracted: 'elevation', 'stage', 'xmomentum' and 'ymomentum'
+def get_mesh_and_quantities_from_file(filename,
+                                      quantities=None,
+                                      verbose=False):
+    """Get and rebuild mesh structure and associated quantities from sww file
+
+    Input:
+        filename - Name os sww file
+        quantities - Names of quantities to load
+
+    Output:
+        mesh - instance of class Interpolate
+               (including mesh and interpolation functionality)
+        quantities - arrays with quantity values at each mesh node
+        time - vector of stored timesteps
+
+    This function is used by e.g.:
+        get_interpolated_quantities_at_polyline_midpoints
+    """
+
+    # FIXME (Ole): Maybe refactor filefunction using this more fundamental code.
+
+    import types
+    from Scientific.IO.NetCDF import NetCDFFile
+    from shallow_water import Domain
+    from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
+
+    if verbose: log.critical('Reading from %s' % filename)
+
+    fid = NetCDFFile(filename, netcdf_mode_r)    # Open existing file for read
+    time = fid.variables['time'][:]    # Time vector
+    time += fid.starttime[0]
+
+    # Get the variables as numeric arrays
+    x = fid.variables['x'][:]                   # x-coordinates of nodes
+    y = fid.variables['y'][:]                   # y-coordinates of nodes
+    elevation = fid.variables['elevation'][:]   # Elevation
+    stage = fid.variables['stage'][:]           # Water level
+    xmomentum = fid.variables['xmomentum'][:]   # Momentum in the x-direction
+    ymomentum = fid.variables['ymomentum'][:]   # Momentum in the y-direction
+
+    # Mesh (nodes (Mx2), triangles (Nx3))
+    nodes = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
+    triangles = fid.variables['volumes'][:]
+
+    # Get geo_reference
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+    except: #AttributeError, e:
+        # Sww files don't have to have a geo_ref
+        geo_reference = None
+
+    if verbose: log.critical('    building mesh from sww file %s' % filename)
+
+    boundary = None
+
+    #FIXME (Peter Row): Should this be in mesh?
+    if fid.smoothing != 'Yes':
+        nodes = nodes.tolist()
+        triangles = triangles.tolist()
+        nodes, triangles, boundary = weed(nodes, triangles, boundary)
+
+    try:
+        mesh = Mesh(nodes, triangles, boundary, geo_reference=geo_reference)
+    except AssertionError, e:
+        fid.close()
+        msg = 'Domain could not be created: %s. "' % e
+        raise DataDomainError, msg
+
+    quantities = {}
+    quantities['elevation'] = elevation
+    quantities['stage'] = stage
+    quantities['xmomentum'] = xmomentum
+    quantities['ymomentum'] = ymomentum
+
+    fid.close()
+
+    return mesh, quantities, time
+
+
+
+##
 # @brief 
 # @parm time 
@@ -1251,6 +1339,4 @@
 
 
-
-
 ##
 # @brief 
@@ -1310,2 +1396,5 @@
 
     return coordinates, volumes, boundary
+
+
+
Index: /trunk/anuga_core/source/anuga/file/test_sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7770)
@@ -8,5 +8,5 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
 from anuga.shallow_water.shallow_water_domain import Domain
-from sww import load_sww_as_domain
+from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file
 
 # boundary functions
@@ -158,4 +158,97 @@
 
 
+    def test_get_mesh_and_quantities_from_sww_file(self):
+        """test_get_mesh_and_quantities_from_sww_file(self):
+        """     
+        
+        # Generate a test sww file with non trivial georeference
+        
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        # Setup
+        from mesh_factory import rectangular
+
+        # Create basic mesh (100m x 5m)
+        width = 5
+        length = 50
+        t_end = 10
+        points, vertices, boundary = rectangular(length, width, 50, 5)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference = Geo_reference(56,308500,6189000))
+
+        domain.set_name('flowtest')
+        swwfile = domain.get_name() + '.sww'
+        domain.set_datadir('.')
+
+        Br = Reflective_boundary(domain)    # Side walls
+        Bd = Dirichlet_boundary([1, 0, 0])  # inflow
+
+        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        for t in domain.evolve(yieldstep=1, finaltime = t_end):
+            pass
+
+        
+        # Read it
+
+        # Get mesh and quantities from sww file
+        X = get_mesh_and_quantities_from_file(swwfile,
+                                              quantities=['elevation',
+                                                          'stage',
+                                                          'xmomentum',
+                                                          'ymomentum'], 
+                                              verbose=False)
+        mesh, quantities, time = X
+        
+
+        # Check that mesh has been recovered
+        assert num.alltrue(mesh.triangles == domain.get_triangles())
+        assert num.allclose(mesh.nodes, domain.get_nodes())
+
+        # Check that time has been recovered
+        assert num.allclose(time, range(t_end+1))
+
+        # Check that quantities have been recovered
+        # (sww files use single precision)
+        z=domain.get_quantity('elevation').get_values(location='unique vertices')
+        assert num.allclose(quantities['elevation'], z)
+
+        for q in ['stage', 'xmomentum', 'ymomentum']:
+            # Get quantity at last timestep
+            q_ref=domain.get_quantity(q).get_values(location='unique vertices')
+
+            #print q,quantities[q]
+            q_sww=quantities[q][-1,:]
+
+            msg = 'Quantity %s failed to be recovered' %q
+            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg
+            
+        # Cleanup
+        os.remove(swwfile)
+        
+        
+
+    def test_weed(self):
+        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
+        volumes1 = [[0,1,2],[3,4,5]]
+        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
+        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
+
+        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
+
+        assert len(points2)==len(coordinates2)
+        for i in range(len(coordinates2)):
+            coordinate = tuple(coordinates2[i])
+            assert points2.has_key(coordinate)
+            points2[coordinate]=i
+
+        for triangle in volumes1:
+            for coordinate in triangle:
+                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
+                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
+
 #################################################################################
 
Index: /trunk/anuga_core/source/anuga/file_conversion/file_conversion.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/file_conversion.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file_conversion/file_conversion.py	(revision 7770)
@@ -488,2 +488,44 @@
 
 
+##
+# @brief 
+# @param filename 
+# @param x 
+# @param y 
+# @param z 
+def write_obj(filename, x, y, z):
+    """Store x,y,z vectors into filename (obj format).
+
+       Vectors are assumed to have dimension (M,3) where
+       M corresponds to the number elements.
+       triangles are assumed to be disconnected
+
+       The three numbers in each vector correspond to three vertices,
+
+       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
+    """
+
+    import os.path
+
+    root, ext = os.path.splitext(filename)
+    if ext == '.obj':
+        FN = filename
+    else:
+        FN = filename + '.obj'
+
+    outfile = open(FN, 'wb')
+    outfile.write("# Triangulation as an obj file\n")
+
+    M, N = x.shape
+    assert N == 3  #Assuming three vertices per element
+
+    for i in range(M):
+        for j in range(N):
+            outfile.write("v %f %f %f\n" % (x[i,j], y[i,j], z[i,j]))
+
+    for i in range(M):
+        base = i * N
+        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
+
+    outfile.close()
+
Index: /trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7770)
@@ -8,9 +8,12 @@
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
 from anuga.shallow_water.sww_file import SWW_file
-from data_manager import extent_sww
+from anuga.file.sww import extent_sww
 from anuga.config import netcdf_float, epsilon, g
 from Scientific.IO.NetCDF import NetCDFFile
 from anuga.file_conversion.file_conversion import tsh2sww, \
                         pmesh_to_domain_instance
+
+from anuga.file.mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
 
 import sys
@@ -1570,5 +1573,5 @@
         #print "sww_file",sww_file
         #dem_file = tempfile.mktemp(".dem")
-        domain = sww2domain(sww_file) ###, dem_file)
+        domain = load_sww_as_domain(sww_file) ###, dem_file)
         domain.check_integrity()
 
Index: /trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py	(revision 7770)
+++ /trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py	(revision 7770)
@@ -0,0 +1,1999 @@
+import numpy as num
+import unittest
+import tempfile
+import os
+from Scientific.IO.NetCDF import NetCDFFile
+
+from anuga.utilities.system_tools import get_pathname_from_package
+from anuga.coordinate_transforms.geo_reference import Geo_reference     
+from anuga.coordinate_transforms.redfearn import redfearn
+from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.file.sts import create_sts_boundary
+from anuga.file.csv_file import load_csv_as_dict, load_csv_as_array
+
+# 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
+
+from urs2sts import urs2sts
+
+# Allow us to use helper methods from this test.
+from anuga.file.test_mux import Test_Mux
+
+class Test_Urs2Sts(Test_Mux):
+    """ A suite of tests to test urs2sts file conversion functions.
+        These tests are quite coarse-grained: converting a file
+        and checking that its headers and some of its contents
+        are correct.
+    """ 
+
+    def test_urs2sts0(self):
+        """
+        Test single source
+        """
+        tide=0
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ha[0]=num.arange(0,time_step_count)
+        ha[1]=num.arange(time_step_count,2*time_step_count)
+        ha[2]=num.arange(2*time_step_count,3*time_step_count)
+        ha[3]=num.arange(3*time_step_count,4*time_step_count)
+        ua=5*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+
+        base_name, files = self.write_mux2(lat_long_points,
+                                      time_step_count, time_step,
+                                      first_tstep, last_tstep,
+                                      depth=gauge_depth,
+                                      ha=ha,
+                                      ua=ua,
+                                      va=va)
+
+        urs2sts(base_name,
+                basename_out=base_name, 
+                mean_stage=tide,verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_name + '.sts'
+
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sts_file)
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        for i in range(4):
+            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1]) 
+            assert num.allclose([x[i],y[i]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        # Set original data used to write mux file to be zero when gauges are
+        #not recdoring
+        ha[0][0]=0.0
+        ha[0][time_step_count-1]=0.0;
+        ha[2][0]=0.0;
+        ua[0][0]=0.0
+        ua[0][time_step_count-1]=0.0;
+        ua[2][0]=0.0;
+        va[0][0]=0.0
+        va[0][time_step_count-1]=0.0;
+        va[2][0]=0.0;
+
+        assert num.allclose(num.transpose(ha),stage)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
+        for i in range(len(lat_long_points)):
+            depth[i]=gauge_depth[i]+tide+ha[i]
+        assert num.allclose(num.transpose(ua*depth),xmomentum) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        assert num.allclose(num.transpose(va*depth),ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert num.allclose(-elevation, gauge_depth)  #Meters
+
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sts_file)
+
+    def test_urs2sts_nonstandard_meridian(self):
+        """
+        Test single source using the meridian from zone 50 as a nonstandard meridian
+        """
+        tide=0
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.,114.5),(-21.,113.5),(-21.,114.), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ha[0]=num.arange(0,time_step_count)
+        ha[1]=num.arange(time_step_count,2*time_step_count)
+        ha[2]=num.arange(2*time_step_count,3*time_step_count)
+        ha[3]=num.arange(3*time_step_count,4*time_step_count)
+        ua=5*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        urs2sts(base_name,
+                basename_out=base_name, 
+                central_meridian=123,
+                mean_stage=tide,
+                verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_name + '.sts'
+
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sts_file)
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        # Check that all coordinate are correctly represented       
+        # Using the non standard projection (50) 
+        for i in range(4):
+            zone, e, n = redfearn(lat_long_points[i][0],
+                                  lat_long_points[i][1],
+                                  central_meridian=123)
+            assert num.allclose([x[i],y[i]], [e,n])
+            assert zone==-1
+        self.delete_mux(files)
+        
+
+
+    def test_urs2sts_individual_sources(self):   
+        """Test that individual sources compare to actual urs output
+           Test that the first recording time is the smallest
+           over waveheight, easting and northing velocity
+        """
+        
+        # Get path where this test is run
+        path = get_pathname_from_package('anuga.shallow_water')        
+        
+        testdir = os.path.join(path, 'urs_test_data')
+        ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt')
+        
+        sources = ['1-z.grd','2-z.grd','3-z.grd']
+        
+        # Start times by source and station taken manually from urs header files
+        time_start_z = num.array([[10.0,11.5,13,14.5,17.7],
+                                  [9.8,11.2,12.7,14.2,17.4],
+                                  [9.5,10.9,12.4,13.9,17.1]])
+
+        time_start_e = time_start_n = time_start_z
+
+        # time step in urs output
+        delta_t = 0.1
+        
+        # Make sts file for each source
+        for k, source_filename in enumerate(sources):
+            source_number = k + 1 # Source numbering starts at 1
+            
+            urs_filenames = os.path.join(testdir, source_filename)
+            weights = [1.]
+            sts_name_out = 'test'
+            
+            urs2sts(urs_filenames,
+                    basename_out=sts_name_out,
+                    ordering_filename=ordering_filename,
+                    weights=weights,
+                    mean_stage=0.0,
+                    verbose=False)
+
+            # Read in sts file for this source file
+            fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r) # Open existing file for read
+            x = fid.variables['x'][:]+fid.xllcorner    # x-coordinates of vertices
+            y = fid.variables['y'][:]+fid.yllcorner    # y-coordinates of vertices
+            elevation = fid.variables['elevation'][:]
+            time=fid.variables['time'][:]+fid.starttime
+
+            # Get quantity data from sts file
+            quantity_names=['stage','xmomentum','ymomentum']
+            quantities = {}
+            for i, name in enumerate(quantity_names):
+                quantities[name] = fid.variables[name][:]
+            
+            # Make sure start time from sts file is the minimum starttime 
+            # across all stations (for this source)
+            #print k, time_start_z[k,:]
+            starttime = min(time_start_z[k, :])
+            sts_starttime = fid.starttime[0]
+            msg = 'sts starttime for source %d was %f. Should have been %f'\
+                %(source_number, sts_starttime, starttime)
+            assert num.allclose(sts_starttime, starttime), msg             
+
+            # For each station, compare urs2sts output to known urs output
+            for j in range(len(x)):
+                index_start_urs = 0
+                index_end_urs = 0
+                index_start = 0
+                index_end = 0
+                count = 0
+
+                # read in urs test data for stage, e and n velocity
+                urs_file_name_z = 'z_'+str(source_number)+'_'+str(j)+'.csv'
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
+                urs_stage = dict['urs_stage']
+                urs_file_name_e = 'e_'+str(source_number)+'_'+str(j)+'.csv'
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
+                urs_e = dict['urs_e']
+                urs_file_name_n = 'n_'+str(source_number)+'_'+str(j)+'.csv'
+                dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
+                urs_n = dict['urs_n']
+
+                # find start and end time for stage             
+                for i in range(len(urs_stage)):
+                    if urs_stage[i] == 0.0:
+                        index_start_urs_z = i+1
+                    if int(urs_stage[i]) == 99 and count <> 1:
+                        count +=1
+                        index_end_urs_z = i
+
+                if count == 0: index_end_urs_z = len(urs_stage)
+
+                start_times_z = time_start_z[source_number-1]
+
+                # start times for easting velocities should be the same as stage
+                start_times_e = time_start_e[source_number-1]
+                index_start_urs_e = index_start_urs_z
+                index_end_urs_e = index_end_urs_z
+
+                # start times for northing velocities should be the same as stage
+                start_times_n = time_start_n[source_number-1]
+                index_start_urs_n = index_start_urs_z
+                index_end_urs_n = index_end_urs_z
+                
+                # Check that actual start time matches header information for stage
+                msg = 'stage start time from urs file is not the same as the '
+                msg += 'header file for source %i and station %i' %(source_number,j)
+                assert num.allclose(index_start_urs_z,start_times_z[j]/delta_t), msg
+
+                msg = 'e velocity start time from urs file is not the same as the '
+                msg += 'header file for source %i and station %i' %(source_number,j)
+                assert num.allclose(index_start_urs_e,start_times_e[j]/delta_t), msg
+
+                msg = 'n velocity start time from urs file is not the same as the '
+                msg += 'header file for source %i and station %i' %(source_number,j)
+                assert num.allclose(index_start_urs_n,start_times_n[j]/delta_t), msg
+                
+                # get index for start and end time for sts quantities
+                index_start_stage = 0
+                index_end_stage = 0
+                count = 0
+                sts_stage = quantities['stage'][:,j]
+                for i in range(len(sts_stage)):
+                    if sts_stage[i] <> 0.0 and count <> 1:
+                        count += 1
+                        index_start_stage = i
+                    if int(sts_stage[i]) == 99 and count <> 1:
+                        count += 1
+                        index_end_stage = i
+
+                index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
+
+                sts_xmom = quantities['xmomentum'][:,j]
+                index_start_x = index_start_stage
+                index_end_x = index_start_x + len(urs_e[index_start_urs_e:index_end_urs_e])
+
+                sts_ymom = quantities['ymomentum'][:,j]
+                index_start_y = index_start_stage
+                index_end_y = index_start_y + len(urs_n[index_start_urs_n:index_end_urs_n])
+
+                # check that urs stage and sts stage are the same
+                msg = 'urs stage is not equal to sts stage for for source %i and station %i' %(source_number,j)
+                assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
+                                    sts_stage[index_start_stage:index_end_stage], 
+                                    rtol=1.0e-6, atol=1.0e-5 ), msg                                
+                                
+                # check that urs e velocity and sts xmomentum are the same
+                msg = 'urs e velocity is not equivalent to sts x momentum for for source %i and station %i' %(source_number,j)
+                assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]),
+                                sts_xmom[index_start_x:index_end_x], 
+                                rtol=1.0e-5, atol=1.0e-4 ), msg
+                
+                # check that urs n velocity and sts ymomentum are the same
+                #print 'urs n velocity', urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j])
+                #print 'sts momentum', sts_ymom[index_start_y:index_end_y]                                                             
+                msg = 'urs n velocity is not equivalent to sts y momentum for source %i and station %i' %(source_number,j)
+                assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
+                                -sts_ymom[index_start_y:index_end_y], 
+                                rtol=1.0e-5, atol=1.0e-4 ), msg
+                                                
+                                
+            fid.close()
+            
+        os.remove(sts_name_out+'.sts')
+
+    def test_urs2sts_combined_sources(self):   
+        """Test that combined sources compare to actual urs output
+           Test that the first recording time is the smallest
+           over waveheight, easting and northing velocity
+        """
+
+        # combined
+        time_start_z = num.array([9.5,10.9,12.4,13.9,17.1])
+        time_start_e = time_start_n = time_start_z
+         
+        # make sts file for combined sources
+        weights = [1., 2., 3.]
+        
+        path = get_pathname_from_package('anuga.shallow_water')        
+                
+        testdir = os.path.join(path, 'urs_test_data')        
+        ordering_filename=os.path.join(testdir, 'thinned_bound_order_test.txt')
+
+        urs_filenames = [os.path.join(testdir,'1-z.grd'),
+                         os.path.join(testdir,'2-z.grd'),
+                         os.path.join(testdir,'3-z.grd')]
+        sts_name_out = 'test'
+        
+        urs2sts(urs_filenames,
+                basename_out=sts_name_out,
+                ordering_filename=ordering_filename,
+                weights=weights,
+                mean_stage=0.0,
+                verbose=False)
+        
+        # read in sts file for combined source
+        fid = NetCDFFile(sts_name_out+'.sts', netcdf_mode_r)    # Open existing file for read
+        x = fid.variables['x'][:]+fid.xllcorner   # x-coordinates of vertices
+        y = fid.variables['y'][:]+fid.yllcorner   # y-coordinates of vertices
+        elevation = fid.variables['elevation'][:]
+        time=fid.variables['time'][:]+fid.starttime
+        
+        
+        # Check that stored permutation is as per default
+        permutation = range(len(x))
+        stored_permutation = fid.variables['permutation'][:]
+        msg = 'Permutation was not stored correctly. I got '
+        msg += str(stored_permutation)
+        assert num.allclose(stored_permutation, permutation), msg        
+
+        # Get quantity data from sts file
+        quantity_names=['stage','xmomentum','ymomentum']
+        quantities = {}
+        for i, name in enumerate(quantity_names):
+            quantities[name] = fid.variables[name][:]
+
+        # For each station, compare urs2sts output to known urs output   
+        delta_t = 0.1
+        
+        # Make sure start time from sts file is the minimum starttime 
+        # across all stations (for this source)
+        starttime = min(time_start_z[:])
+        sts_starttime = fid.starttime[0]
+        msg = 'sts starttime was %f. Should have been %f'\
+            %(sts_starttime, starttime)
+        assert num.allclose(sts_starttime, starttime), msg
+    
+        #stations = [1,2,3]
+        #for j in stations: 
+        for j in range(len(x)):
+            index_start_urs_z = 0
+            index_end_urs_z = 0
+            index_start_urs_e = 0
+            index_end_urs_e = 0
+            index_start_urs_n = 0
+            index_end_urs_n = 0
+            count = 0
+
+            # read in urs test data for stage, e and n velocity
+            urs_file_name_z = 'z_combined_'+str(j)+'.csv'
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_z))
+            urs_stage = dict['urs_stage']
+            urs_file_name_e = 'e_combined_'+str(j)+'.csv'
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_e))
+            urs_e = dict['urs_e']
+            urs_file_name_n = 'n_combined_'+str(j)+'.csv'
+            dict = load_csv_as_array(os.path.join(testdir, urs_file_name_n))
+            urs_n = dict['urs_n']
+
+            # find start and end time for stage         
+            for i in range(len(urs_stage)):
+                if urs_stage[i] == 0.0:
+                    index_start_urs_z = i+1
+                if int(urs_stage[i]) == 99 and count <> 1:
+                    count +=1
+                    index_end_urs_z = i
+
+            if count == 0: index_end_urs_z = len(urs_stage)
+
+            start_times_z = time_start_z[j]
+
+            start_times_e = time_start_e[j]
+            index_start_urs_e = index_start_urs_z
+
+            start_times_n = time_start_n[j]
+            index_start_urs_n = index_start_urs_z
+               
+            # Check that actual start time matches header information for stage
+            msg = 'stage start time from urs file is not the same as the '
+            msg += 'header file at station %i' %(j)
+            assert num.allclose(index_start_urs_z,start_times_z/delta_t), msg
+
+            msg = 'e velocity start time from urs file is not the same as the '
+            msg += 'header file at station %i' %(j)
+            assert num.allclose(index_start_urs_e,start_times_e/delta_t), msg
+
+            msg = 'n velocity start time from urs file is not the same as the '
+            msg += 'header file at station %i' %(j)
+            assert num.allclose(index_start_urs_n,start_times_n/delta_t), msg
+                
+            # get index for start and end time for sts quantities
+            index_start_stage = 0
+            index_end_stage = 0
+            index_start_x = 0
+            index_end_x = 0
+            index_start_y = 0
+            index_end_y = 0
+            count = 0
+            count1 = 0
+            sts_stage = quantities['stage'][:,j]
+            for i in range(len(sts_stage)):
+                if sts_stage[i] <> 0.0 and count <> 1:
+                    count += 1
+                    index_start_stage = i
+                if int(urs_stage[i]) == 99 and count <> 1:
+                    count +=1
+                    index_end_stage = i
+                
+            index_end_stage = index_start_stage + len(urs_stage[index_start_urs_z:index_end_urs_z])
+
+            index_start_x = index_start_stage
+            index_end_x = index_start_x + len(urs_stage[index_start_urs_e:index_end_urs_e])
+            sts_xmom = quantities['ymomentum'][:,j]
+
+            index_start_y = index_start_stage
+            index_end_y = index_start_y + len(urs_stage[index_start_urs_n:index_end_urs_n])
+            sts_ymom = quantities['ymomentum'][:,j]
+
+            # check that urs stage and sts stage are the same
+            msg = 'urs stage is not equal to sts stage for station %i' %j
+            #print 'urs stage', urs_stage[index_start_urs_z:index_end_urs_z]
+            #print 'sts stage', sts_stage[index_start_stage:index_end_stage]
+            #print 'diff', max(urs_stage[index_start_urs_z:index_end_urs_z]-sts_stage[index_start_stage:index_end_stage])
+            #print 'index', index_start_stage, index_end_stage, len(sts_stage)
+            assert num.allclose(urs_stage[index_start_urs_z:index_end_urs_z],
+                            sts_stage[index_start_stage:index_end_stage], 
+                                rtol=1.0e-5, atol=1.0e-4 ), msg                                
+                                
+            # check that urs e velocity and sts xmomentum are the same          
+            msg = 'urs e velocity is not equivalent to sts xmomentum for station %i' %j
+            assert num.allclose(urs_e[index_start_urs_e:index_end_urs_e]*(urs_stage[index_start_urs_e:index_end_urs_e]-elevation[j]),
+                            sts_xmom[index_start_x:index_end_x], 
+                            rtol=1.0e-5, atol=1.0e-4 ), msg
+                
+            # check that urs n velocity and sts ymomentum are the same                            
+            msg = 'urs n velocity is not equivalent to sts ymomentum for station %i' %j
+            assert num.allclose(urs_n[index_start_urs_n:index_end_urs_n]*(urs_stage[index_start_urs_n:index_end_urs_n]-elevation[j]),
+                            sts_ymom[index_start_y:index_end_y], 
+                            rtol=1.0e-5, atol=1.0e-4 ), msg
+
+        fid.close()
+        
+        os.remove(sts_name_out+'.sts')
+        
+        
+        
+    def test_urs2sts_ordering(self):
+        """Test multiple sources with ordering file
+        """
+        
+        tide = 0.35
+        time_step_count = 6 # I made this a little longer (Ole)
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ha[0]=num.arange(0,time_step_count)
+        ha[1]=num.arange(time_step_count,2*time_step_count)
+        ha[2]=num.arange(2*time_step_count,3*time_step_count)
+        ha[3]=num.arange(3*time_step_count,4*time_step_count)
+        ua=5*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+
+        # Create two identical mux files to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha,
+                                             ua=ua,
+                                             va=va)
+
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha,
+                                               ua=ua,
+                                               va=va)
+
+                                               
+        # Create ordering file
+        permutation = [3,0,2]
+
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][1], 
+                                             lat_long_points[index][0]))
+        order_fid.close()
+        
+            
+
+                                               
+        # Call urs2sts with multiple mux files
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                weights=[1.0, 1.0],
+                mean_stage=tide,
+                verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_nameI + '.sts'
+
+        #Let's interrogate the sts file
+        # Note, the sts info is not gridded.  It is point data.
+        fid = NetCDFFile(sts_file)
+        
+        # Check that original indices have been stored
+        stored_permutation = fid.variables['permutation'][:]
+        msg = 'Permutation was not stored correctly. I got '
+        msg += str(stored_permutation)
+        assert num.allclose(stored_permutation, permutation), msg
+        
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        #print
+        #print x
+        #print y
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+
+            #print i, [x[i],y[i]], [e,n]
+            assert num.allclose([x[i],y[i]], [e,n])
+            
+                        
+        # Check the time vector
+        times = fid.variables['time'][:]
+
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+                        
+
+        # Check sts values
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        # Set original data used to write mux file to be zero when gauges are
+        # not recdoring
+        
+        ha[0][0]=0.0
+        ha[0][time_step_count-1]=0.0
+        ha[2][0]=0.0
+        ua[0][0]=0.0
+        ua[0][time_step_count-1]=0.0
+        ua[2][0]=0.0
+        va[0][0]=0.0
+        va[0][time_step_count-1]=0.0
+        va[2][0]=0.0;
+
+        
+        # The stage stored in the .sts file should be the sum of the stage
+        # in the two mux2 files because both have weights = 1. In this case
+        # the mux2 files are the same so stage == 2.0 * ha
+        #print 2.0*num.transpose(ha) - stage 
+        
+        ha_permutation = num.take(ha, permutation, axis=0) 
+        ua_permutation = num.take(ua, permutation, axis=0)
+        va_permutation = num.take(va, permutation, axis=0)
+        gauge_depth_permutation = num.take(gauge_depth, permutation, axis=0)
+
+        
+        assert num.allclose(2.0*num.transpose(ha_permutation)+tide, stage)  # Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
+        for i in range(len(lat_long_points)):
+            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
+            #2.0*ha necessary because using two files with weights=1 are used
+            
+        depth_permutation = num.take(depth, permutation, axis=0)                     
+        
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert num.allclose(2.0*num.transpose(ua_permutation*depth_permutation), xmomentum) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        assert num.allclose(2.0*num.transpose(va_permutation*depth_permutation), ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert num.allclose(-gauge_depth_permutation, elevation)  #Meters
+
+        fid.close()
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
+        os.remove(sts_file)
+        
+
+        
+        
+        
+    def Xtest_urs2sts_ordering_exception(self):
+        """Test that inconsistent lats and lons in ordering file are caught.
+        """
+        
+        tide=0
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ha[0]=num.arange(0,time_step_count)
+        ha[1]=num.arange(time_step_count,2*time_step_count)
+        ha[2]=num.arange(2*time_step_count,3*time_step_count)
+        ha[3]=num.arange(3*time_step_count,4*time_step_count)
+        ua=5*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+
+        # Create two identical mux files to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha,
+                                             ua=ua,
+                                             va=va)
+
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha,
+                                               ua=ua,
+                                               va=va)
+
+                                               
+        # Create ordering file
+        permutation = [3,0,2]
+
+        # Do it wrongly and check that exception is being raised
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][0], 
+                                             lat_long_points[index][1]))
+        order_fid.close()
+        
+        try:
+            urs2sts([base_nameI, base_nameII], 
+                    basename_out=base_nameI, 
+                    ordering_filename=ordering_filename,
+                    weights=[1.0, 1.0],
+                    mean_stage=tide,
+                    verbose=False)  
+            os.remove(ordering_filename)            
+        except:
+            pass
+        else:
+            msg = 'Should have caught wrong lat longs'
+            raise Exception, msg
+
+        
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
+
+        
+
+        
+    def test_urs2sts_ordering_different_sources(self):
+        """Test multiple sources with ordering file, different source files and weights.
+           This test also has more variable values than the previous ones
+        """
+        
+        tide = 1.5
+        time_step_count = 10
+        time_step = 0.2
+        
+        times_ref = num.arange(0, time_step_count*time_step, time_step)
+        #print 'time vector', times_ref
+        
+        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
+        n = len(lat_long_points)
+        
+        # Create non-trivial weights
+        #weights = [0.8, 1.5] # OK
+        #weights = [0.8, 10.5] # Fail (up to allclose tolerance)
+        #weights = [10.5, 10.5] # OK
+        #weights = [0.0, 10.5] # OK
+        #weights = [0.8, 0.] # OK                
+        #weights = [8, 0.1] # OK                        
+        #weights = [0.8, 10.0] # OK                                
+        #weights = [0.8, 10.6] # OK           
+        weights = [3.8, 7.6] # OK                   
+        #weights = [0.5, 0.5] # OK                           
+        #weights = [2., 2.] # OK                            
+        #weights = [0.0, 0.5] # OK                                          
+        #weights = [1.0, 1.0] # OK                                                  
+        
+        
+        # Create different timeseries starting and ending at different times 
+        first_tstep=num.ones(n,num.int)
+        first_tstep[0]+=2   # Point 0 starts at 2
+        first_tstep[1]+=4   # Point 1 starts at 4        
+        first_tstep[2]+=3   # Point 2 starts at 3
+        
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        last_tstep[0]-=1    # Point 0 ends 1 step early
+        last_tstep[1]-=2    # Point 1 ends 2 steps early                
+        last_tstep[4]-=3    # Point 4 ends 3 steps early        
+        
+        #print
+        #print 'time_step_count', time_step_count
+        #print 'time_step', time_step
+        #print 'first_tstep', first_tstep
+        #print 'last_tstep', last_tstep                
+        
+        
+        # Create varying elevation data (positive values for seafloor)
+        gauge_depth=20*num.ones(n,num.float)
+        for i in range(n):
+            gauge_depth[i] += i**2
+            
+        #print 'gauge_depth', gauge_depth
+        
+        # Create data to be written to first mux file        
+        ha0=2*num.ones((n,time_step_count),num.float)
+        ha0[0]=num.arange(0,time_step_count)
+        ha0[1]=num.arange(time_step_count,2*time_step_count)
+        ha0[2]=num.arange(2*time_step_count,3*time_step_count)
+        ha0[3]=num.arange(3*time_step_count,4*time_step_count)
+        ua0=5*num.ones((n,time_step_count),num.float)
+        va0=-10*num.ones((n,time_step_count),num.float)
+
+        # Ensure data used to write mux file to be zero when gauges are
+        # not recording
+        for i in range(n):
+             # For each point
+             
+             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
+                 # For timesteps before and after recording range
+                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                  
+
+
+                 
+        #print 
+        #print 'using varying start and end time'
+        #print 'ha0', ha0[4]
+        #print 'ua0', ua0
+        #print 'va0', va0        
+        
+        # Write first mux file to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha0,
+                                             ua=ua0,
+                                             va=va0)
+
+                                             
+                                             
+        # Create data to be written to second mux file        
+        ha1=num.ones((n,time_step_count),num.float)
+        ha1[0]=num.sin(times_ref)
+        ha1[1]=2*num.sin(times_ref - 3)
+        ha1[2]=5*num.sin(4*times_ref)
+        ha1[3]=num.sin(times_ref)
+        ha1[4]=num.sin(2*times_ref-0.7)
+                
+        ua1=num.zeros((n,time_step_count),num.float)
+        ua1[0]=3*num.cos(times_ref)        
+        ua1[1]=2*num.sin(times_ref-0.7)   
+        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
+        ua1[4]=2*num.ones(time_step_count)
+        
+        va1=num.zeros((n,time_step_count),num.float)
+        va1[0]=2*num.cos(times_ref-0.87)        
+        va1[1]=3*num.ones(time_step_count, num.int)       #array default#
+        va1[3]=2*num.sin(times_ref-0.71)        
+        
+        
+        # Ensure data used to write mux file to be zero when gauges are
+        # not recording
+        for i in range(n):
+             # For each point
+             
+             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
+                 # For timesteps before and after recording range
+                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                  
+
+
+        #print 
+        #print 'using varying start and end time'
+        #print 'ha1', ha1[4]
+        #print 'ua1', ua1
+        #print 'va1', va1        
+                                             
+                                             
+        # Write second mux file to be combined by urs2sts                                             
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha1,
+                                               ua=ua1,
+                                               va=va1)
+
+                                               
+        # Create ordering file
+        permutation = [4,0,2]
+
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][1], 
+                                             lat_long_points[index][0]))
+        order_fid.close()
+        
+            
+
+        #------------------------------------------------------------
+        # Now read the mux files one by one without weights and test
+        
+        # Call urs2sts with mux file #0
+        urs2sts([base_nameI], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                mean_stage=tide,
+                verbose=False)
+
+        # Now read the sts file and check that values have been stored correctly.
+        sts_file = base_nameI + '.sts'
+        fid = NetCDFFile(sts_file)
+        
+
+        # Check that original indices have been stored
+        stored_permutation = fid.variables['permutation'][:]
+        msg = 'Permutation was not stored correctly. I got '
+        msg += str(stored_permutation)
+        assert num.allclose(stored_permutation, permutation), msg
+        
+
+        
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+
+            #print i, [x[i],y[i]], [e,n]
+            assert num.allclose([x[i],y[i]], [e,n])
+            
+                        
+        # Check the time vector
+        times = fid.variables['time'][:]
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_ref))
+                        
+
+        # Check sts values for mux #0
+        stage0 = fid.variables['stage'][:].copy()
+        xmomentum0 = fid.variables['xmomentum'][:].copy()
+        ymomentum0 = fid.variables['ymomentum'][:].copy()
+        elevation0 = fid.variables['elevation'][:].copy()
+
+        
+        #print 'stage', stage0
+        #print 'xmomentum', xmomentum0
+        #print 'ymomentum', ymomentum0        
+        #print 'elevation', elevation0
+        
+        # The quantities stored in the .sts file should be the weighted sum of the 
+        # quantities written to the mux2 files subject to the permutation vector.
+        
+        ha_ref = num.take(ha0, permutation, axis=0)
+        ua_ref = num.take(ua0, permutation, axis=0)        
+        va_ref = num.take(va0, permutation, axis=0)                
+
+        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                      
+        
+        assert num.allclose(num.transpose(ha_ref)+tide, stage0)  # Meters
+        
+        
+        
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
+        for i in range(len(permutation)):
+            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
+
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum0) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        
+        
+        #print transpose(va_ref*depth_ref)
+        #print ymomentum
+        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum0)        
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert num.allclose(-gauge_depth_ref, elevation0) 
+
+        fid.close()
+        os.remove(sts_file)
+        
+        
+
+        
+        # Call urs2sts with mux file #1
+        urs2sts([base_nameII], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                mean_stage=tide,
+                verbose=False)
+
+        # Now read the sts file and check that values have been stored correctly.
+        sts_file = base_nameI + '.sts'
+        fid = NetCDFFile(sts_file)
+        
+        
+        # Check that original indices have been stored
+        stored_permutation = fid.variables['permutation'][:]
+        msg = 'Permutation was not stored correctly. I got '
+        msg += str(stored_permutation)
+        assert num.allclose(stored_permutation, permutation), msg
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+
+            #print i, [x[i],y[i]], [e,n]
+            assert num.allclose([x[i],y[i]], [e,n])
+            
+                        
+        # Check the time vector
+        times = fid.variables['time'][:]
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_ref))
+                        
+
+        # Check sts values for mux #1 
+        stage1 = fid.variables['stage'][:].copy()
+        xmomentum1 = fid.variables['xmomentum'][:].copy()
+        ymomentum1 = fid.variables['ymomentum'][:].copy()
+        elevation1 = fid.variables['elevation'][:].copy()
+
+        
+        #print 'stage', stage1
+        #print 'xmomentum', xmomentum1
+        #print 'ymomentum', ymomentum1       
+        #print 'elevation', elevation1
+        
+        # The quantities stored in the .sts file should be the weighted sum of the 
+        # quantities written to the mux2 files subject to the permutation vector.
+        
+        ha_ref = num.take(ha1, permutation, axis=0)
+        ua_ref = num.take(ua1, permutation, axis=0)        
+        va_ref = num.take(va1, permutation, axis=0)                
+
+        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
+
+
+        #print 
+        #print stage1
+        #print transpose(ha_ref)+tide - stage1
+        
+
+        assert num.allclose(num.transpose(ha_ref)+tide, stage1)  # Meters
+        #import sys; sys.exit()
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
+        for i in range(len(permutation)):
+            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
+
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum1) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        
+        
+        #print transpose(va_ref*depth_ref)
+        #print ymomentum
+        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum1)        
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert num.allclose(-gauge_depth_ref, elevation1) 
+
+        fid.close()
+        os.remove(sts_file)
+        
+        #----------------------
+        # Then read the mux files together and test
+        
+                                               
+        # Call urs2sts with multiple mux files
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                weights=weights,
+                mean_stage=tide,
+                verbose=False)
+
+        # Now read the sts file and check that values have been stored correctly.
+        sts_file = base_nameI + '.sts'
+        fid = NetCDFFile(sts_file)
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+
+            #print i, [x[i],y[i]], [e,n]
+            assert num.allclose([x[i],y[i]], [e,n])
+            
+                        
+        # Check the time vector
+        times = fid.variables['time'][:]
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_ref))
+                        
+
+        # Check sts values
+        stage = fid.variables['stage'][:].copy()
+        xmomentum = fid.variables['xmomentum'][:].copy()
+        ymomentum = fid.variables['ymomentum'][:].copy()
+        elevation = fid.variables['elevation'][:].copy()
+
+        
+        #print 'stage', stage
+        #print 'elevation', elevation
+        
+        # The quantities stored in the .sts file should be the weighted sum of the 
+        # quantities written to the mux2 files subject to the permutation vector.
+        
+        ha_ref = (weights[0]*num.take(ha0, permutation, axis=0)
+                  + weights[1]*num.take(ha1, permutation, axis=0))
+        ua_ref = (weights[0]*num.take(ua0, permutation, axis=0)
+                  + weights[1]*num.take(ua1, permutation, axis=0))
+        va_ref = (weights[0]*num.take(va0, permutation, axis=0)
+                  + weights[1]*num.take(va1, permutation, axis=0))
+
+        gauge_depth_ref = num.take(gauge_depth, permutation, axis=0)                         
+
+
+        #print 
+        #print stage
+        #print transpose(ha_ref)+tide - stage
+
+        assert num.allclose(num.transpose(ha_ref)+tide, stage)  # Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth_ref = num.zeros((len(permutation), time_step_count), num.float)
+        for i in range(len(permutation)):
+            depth_ref[i]=gauge_depth_ref[i]+tide+ha_ref[i]
+
+            
+        
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert num.allclose(num.transpose(ua_ref*depth_ref), xmomentum) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        
+        
+        #print transpose(va_ref*depth_ref)
+        #print ymomentum
+
+        assert num.allclose(num.transpose(va_ref*depth_ref), ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert num.allclose(-gauge_depth_ref, elevation)  #Meters
+
+        fid.close()
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
+        os.remove(sts_file)
+        
+        #---------------
+        # "Manually" add the timeseries up with weights and test
+        # Tide is discounted from individual results and added back in       
+        #
+
+        stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
+        assert num.allclose(stage_man, stage)
+                
+        
+    def test_file_boundary_stsI(self):
+        """test_file_boundary_stsI(self):
+        """
+        
+        # FIXME (Ole): These tests should really move to 
+        # test_generic_boundaries.py
+        
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide = 0.37
+        time_step_count = 5
+        time_step = 2
+        lat_long_points =bounding_polygon[0:3]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+
+        h = 20        
+        w = 2
+        u = 10
+        v = -10
+        gauge_depth=h*num.ones(n,num.float)
+        ha=w*num.ones((n,time_step_count),num.float)
+        ua=u*num.ones((n,time_step_count),num.float)
+        va=v*num.ones((n,time_step_count),num.float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1])
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=bounding_polygon)
+        Br = Reflective_boundary(domain_fbound)
+        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])        
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        
+        # Check boundary object evaluates as it should
+        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
+            if B is Bf:
+            
+                qf = B.evaluate(vol_id, edge_id)  # File boundary
+                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
+
+                assert num.allclose(qf, qd) 
+                
+        
+        # Evolve
+        finaltime=time_step*(time_step_count-1)
+        yieldstep=time_step
+        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
+
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+                                                   
+            D = domain_fbound
+            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
+
+            # Check that file boundary object has populated 
+            # boundary array correctly  
+            # FIXME (Ole): Do this for the other tests too!
+            for j, val in enumerate(D.get_quantity('stage').boundary_values):
+            
+                (vol_id, edge_id), B = D.boundary_objects[j]
+                if isinstance(B, File_boundary):
+                    #print j, val
+                    assert num.allclose(val, w + tide)
+
+
+        
+        domain_drchlt = Domain(meshname)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
+
+        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+
+        #print domain_fbound.quantities['stage'].vertex_values
+        #print domain_drchlt.quantities['stage'].vertex_values
+                    
+        assert num.allclose(temp_fbound,temp_drchlt)
+        
+        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
+                            domain_drchlt.quantities['stage'].vertex_values)
+                        
+        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
+                            domain_drchlt.quantities['xmomentum'].vertex_values) 
+                        
+        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
+                            domain_drchlt.quantities['ymomentum'].vertex_values)
+        
+        
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
+                
+        
+    def test_file_boundary_stsI_beyond_model_time(self):
+        """test_file_boundary_stsI(self):
+        
+        Test that file_boundary can work when model time
+        exceeds available data.
+        This is optional and requires the user to specify a default 
+        boundary object.
+        """
+        
+        # Don't do warnings in unit test
+        import warnings
+        warnings.simplefilter('ignore')
+        
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
+                          [6.02,97.02],[6.00,97.02]]
+        tide = 0.37
+        time_step_count = 5
+        time_step = 2
+        lat_long_points = bounding_polygon[0:3]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+
+        h = 20        
+        w = 2
+        u = 10
+        v = -10
+        gauge_depth=h*num.ones(n,num.float)
+        ha=w*num.ones((n,time_step_count),num.float)
+        ua=u*num.ones((n,time_step_count),num.float)
+        va=v*num.ones((n,time_step_count),num.float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,\
+            bounding_polygon[i][0],\
+            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
+                                            bounding_polygon[i][1])
+                                            
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        Br = Reflective_boundary(domain_fbound)
+        Bd = Dirichlet_boundary([w+tide, u*(w+h+tide), v*(w+h+tide)])        
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=bounding_polygon,
+                           default_boundary=Bd) # Condition to be used 
+                                                # if model time exceeds 
+                                                # available data
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        
+        # Check boundary object evaluates as it should
+        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
+            if B is Bf:
+            
+                qf = B.evaluate(vol_id, edge_id)  # File boundary
+                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
+
+                assert num.allclose(qf, qd) 
+                
+        
+        # Evolve
+        data_finaltime = time_step*(time_step_count-1)
+        finaltime = data_finaltime + 10 # Let model time exceed available data
+        yieldstep = time_step
+        temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float)
+
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+                                                   
+            D = domain_fbound
+            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
+
+            # Check that file boundary object has populated 
+            # boundary array correctly  
+            # FIXME (Ole): Do this for the other tests too!
+            for j, val in enumerate(D.get_quantity('stage').boundary_values):
+            
+                (vol_id, edge_id), B = D.boundary_objects[j]
+                if isinstance(B, File_boundary):
+                    #print j, val
+                    assert num.allclose(val, w + tide)
+
+
+        domain_drchlt = Domain(meshname)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
+
+        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+
+        #print domain_fbound.quantities['stage'].vertex_values
+        #print domain_drchlt.quantities['stage'].vertex_values
+                    
+        assert num.allclose(temp_fbound,temp_drchlt)
+        
+        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
+                            domain_drchlt.quantities['stage'].vertex_values)
+                        
+        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
+                            domain_drchlt.quantities['xmomentum'].vertex_values) 
+                        
+        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
+                            domain_drchlt.quantities['ymomentum'].vertex_values) 
+        
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
+                
+        
+    def test_field_boundary_stsI_beyond_model_time(self):
+        """test_field_boundary(self):
+        
+        Test that field_boundary can work when model time
+        exceeds available data whilst adjusting mean_stage.
+        
+        """
+        
+        # Don't do warnings in unit test
+        import warnings
+        warnings.simplefilter('ignore')
+        
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
+                          [6.02,97.02],[6.00,97.02]]
+        tide = 0.37
+        time_step_count = 5
+        time_step = 2
+        lat_long_points = bounding_polygon[0:3]
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+
+        h = 20        
+        w = 2
+        u = 10
+        v = -10
+        gauge_depth=h*num.ones(n,num.float)
+        ha=w*num.ones((n,time_step_count),num.float)
+        ua=u*num.ones((n,time_step_count),num.float)
+        va=v*num.ones((n,time_step_count),num.float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=0.0, # Deliberately let Field_boundary do the adjustment
+                verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,\
+            bounding_polygon[i][0],\
+            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
+                                            bounding_polygon[i][1])
+                                            
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        Br = Reflective_boundary(domain_fbound)
+        Bd = Dirichlet_boundary([w+tide, u*(w+h), v*(w+h)])
+        Bdefault = Dirichlet_boundary([w, u*(w+h), v*(w+h)])        
+                
+        Bf = Field_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           mean_stage=tide, # Field boundary to adjust for tide
+                           boundary_polygon=bounding_polygon,
+                           default_boundary=Bdefault) # Condition to be used 
+                                                      # if model time exceeds 
+                                                      # available data
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        
+        # Check boundary object evaluates as it should
+        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
+            if B is Bf:
+            
+                qf = B.evaluate(vol_id, edge_id)  # Field boundary
+                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
+                
+                msg = 'Got %s, should have been %s' %(qf, qd)
+                assert num.allclose(qf, qd), msg 
+                
+        # Evolve
+        data_finaltime = time_step*(time_step_count-1)
+        finaltime = data_finaltime + 10 # Let model time exceed available data
+        yieldstep = time_step
+        temp_fbound=num.zeros(int(finaltime/yieldstep)+1, num.float)
+
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+                                                   
+            D = domain_fbound
+            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
+
+            # Check that file boundary object has populated 
+            # boundary array correctly  
+            # FIXME (Ole): Do this for the other tests too!
+            for j, val in enumerate(D.get_quantity('stage').boundary_values):
+            
+                (vol_id, edge_id), B = D.boundary_objects[j]
+                if isinstance(B, Field_boundary):
+                    msg = 'Got %f should have been %f' %(val, w+tide)
+                    assert num.allclose(val, w + tide), msg
+
+
+    def test_file_boundary_stsII(self):
+        """test_file_boundary_stsII(self):
+        
+         mux2 file has points not included in boundary
+         mux2 gauges are not stored with the same order as they are 
+         found in bounding_polygon. This does not matter as long as bounding
+         polygon passed to file_function contains the mux2 points needed (in
+         the correct order).
+         """
+         
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]]
+        tide = -2.20 
+        time_step_count = 20
+        time_step = 2
+        lat_long_points=bounding_polygon[0:2]
+        lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1])
+        lat_long_points.insert(0,[6.0,97.01])
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ua=10*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count,
+                                           time_step,
+                                           first_tstep,
+                                           last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,\
+            bounding_polygon[i][0],\
+            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
+                                            bounding_polygon[i][1])
+            
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        # have to change boundary tags from last example because now bounding
+        # polygon starts in different place.
+        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
+                         maximum_triangle_area=extent_res,filename=meshname,
+                         interior_regions=interior_regions,verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        Bf = File_boundary(sts_file+'.sts',
+                           domain_fbound,
+                           boundary_polygon=bounding_polygon)
+        Br = Reflective_boundary(domain_fbound)
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        finaltime=time_step*(time_step_count-1)
+        yieldstep=time_step
+        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
+        
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step = False)):
+            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
+        
+        domain_drchlt = Domain(meshname)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
+
+        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step = False)):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+
+
+        assert num.allclose(temp_fbound,temp_drchlt)            
+            
+        #print domain_fbound.quantities['stage'].vertex_values
+        #print domain_drchlt.quantities['stage'].vertex_values
+                    
+            
+        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
+                            domain_drchlt.quantities['stage'].vertex_values)
+                        
+        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
+                            domain_drchlt.quantities['xmomentum'].vertex_values)
+                        
+        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
+                            domain_drchlt.quantities['ymomentum'].vertex_values)
+            
+            
+
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
+
+        
+        
+    def test_file_boundary_stsIII_ordering(self):
+        """test_file_boundary_stsIII_ordering(self):
+        Read correct points from ordering file and apply sts to boundary
+        """
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
+                          [6.02,97.02],[6.00,97.02]]
+        tide = 3.0
+        time_step_count = 50
+        time_step = 2
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        gauge_depth=20*num.ones(n,num.float)
+        ha=2*num.ones((n,time_step_count),num.float)
+        ua=10*num.ones((n,time_step_count),num.float)
+        va=-10*num.ones((n,time_step_count),num.float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count,
+                                           time_step,
+                                           first_tstep,
+                                           last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        # Write order file
+        file_handle, order_base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(order_base_name)
+        d=","
+        order_file=order_base_name+'order.txt'
+        fid=open(order_file,'w')
+        
+        # Write Header
+        header='index, longitude, latitude\n'
+        fid.write(header)
+        indices=[3,0,1]
+        for i in indices:
+            line=str(i)+d+str(lat_long_points[i][1])+d+\
+                str(lat_long_points[i][0])+"\n"
+            fid.write(line)
+        fid.close()
+
+        sts_file=base_name
+        urs2sts(base_name,
+                basename_out=sts_file,
+                ordering_filename=order_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+
+        assert(os.access(base_name, os.R_OK))
+
+        boundary_polygon = create_sts_boundary(base_name)
+
+        os.remove(order_file)
+
+        # Append the remaining part of the boundary polygon to be defined by
+        # the user
+        bounding_polygon_utm=[]
+        for point in bounding_polygon:
+            zone,easting,northing=redfearn(point[0],point[1])
+            bounding_polygon_utm.append([easting,northing])
+
+        boundary_polygon.append(bounding_polygon_utm[3])
+        boundary_polygon.append(bounding_polygon_utm[4])
+
+        plot=False
+        if plot:
+            from pylab import plot,show,axis
+            boundary_polygon=ensure_numeric(boundary_polygon)
+            bounding_polygon_utm=ensure_numeric(bounding_polygon_utm)
+            #lat_long_points=ensure_numeric(lat_long_points)
+            #plot(lat_long_points[:,0],lat_long_points[:,1],'o')
+            plot(boundary_polygon[:,0], boundary_polygon[:,1],'d')
+            plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1],'o')
+            show()
+
+        assert num.allclose(bounding_polygon_utm,boundary_polygon)
+
+
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        
+        # have to change boundary tags from last example because now bounding
+        # polygon starts in different place.
+        create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags,
+                         maximum_triangle_area=extent_res,filename=meshname,
+                         interior_regions=interior_regions,verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        Bf = File_boundary(sts_file+'.sts',
+                           domain_fbound,
+                           boundary_polygon=boundary_polygon)
+        Br = Reflective_boundary(domain_fbound)
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        finaltime=time_step*(time_step_count-1)
+        yieldstep=time_step
+        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
+    
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step = False)):
+            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
+    
+        
+        domain_drchlt = Domain(meshname)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+        Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
+        
+        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step = False)):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+
+        
+        #print domain_fbound.quantities['stage'].vertex_values
+        #print domain_drchlt.quantities['stage'].vertex_values
+                    
+        assert num.allclose(temp_fbound,temp_drchlt)
+
+        
+        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
+                            domain_drchlt.quantities['stage'].vertex_values)
+                        
+        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
+                            domain_drchlt.quantities['xmomentum'].vertex_values)                        
+                        
+        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
+                            domain_drchlt.quantities['ymomentum'].vertex_values)
+        
+        # Use known Dirichlet condition (if sufficient timesteps have been taken)
+
+        #FIXME: What do these assertions test? Also do they assume tide =0
+        #print domain_fbound.quantities['stage'].vertex_values
+        #assert allclose(domain_drchlt.quantities['stage'].vertex_values[6], 2)        
+        #assert allclose(domain_fbound.quantities['stage'].vertex_values[6], 2)
+        
+        
+
+        try:
+            os.remove(sts_file+'.sts')
+        except:
+            # Windoze can't remove this file for some reason 
+            pass
+        
+        os.remove(meshname)
+        
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Urs2Sts,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/file_conversion/test_urs2sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/test_urs2sww.py	(revision 7770)
+++ /trunk/anuga_core/source/anuga/file_conversion/test_urs2sww.py	(revision 7770)
@@ -0,0 +1,674 @@
+
+# External modules
+from Scientific.IO.NetCDF import NetCDFFile
+import sys
+import unittest
+import numpy as num
+import copy
+import os
+
+# ANUGA modules
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     write_NetCDF_georeference
+from anuga.coordinate_transforms.redfearn import redfearn     
+
+from urs2sww import urs2sww
+import urs
+
+from anuga.file.mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
+
+from anuga.geospatial_data.geospatial_data import ensure_absolute
+from anuga.utilities.numerical_tools import ensure_numeric  
+
+# use helper methods from other unit test
+from anuga.file.test_mux import Test_Mux
+
+
+
+class Test_Dem2Pts(Test_Mux):
+    """ A suite of tests to test urs2sww file conversion functions.
+        These tests are quite coarse-grained: converting a file
+        and checking that its headers and some of its contents
+        are correct.
+    """ 
+
+    def setUp(self):
+        self.verbose = False # change this to output more debug info
+
+  
+    def test_urs_ungridded2swwIII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 2
+        time_step = 400
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide, origin =(50,23432,4343))
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+        # = n*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_x = n*(e+tide+n)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert num.allclose(answer_x, actual_x)  #Meters
+        
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+        # = e*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_y = -1*e*(e+tide+n)
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert num.allclose(answer_y, actual_y)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+
+        
+    def test_urs_ungridded_hole (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-20.5, 114.5],
+                    [-20.6, 114.6],
+                    [-20.5, 115.],
+                    [-20.6, 115.],
+                    [-20.5, 115.5],
+                    [-20.6, 115.4],
+                    
+                    [-21., 114.5],
+                    [-21., 114.6],
+                    [-21., 115.5],
+                    [-21., 115.4],
+                    
+                    [-21.5, 114.5],
+                    [-21.4, 114.6],
+                    [-21.5, 115.],
+                    [-21.4, 115.],
+                    [-21.5, 115.5],
+                    [-21.4, 115.4]
+                    ]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        #Easting:  292110.784  Northing: 7676551.710 
+        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 '' 
+
+        urs2sww(base_name, mean_stage=-240992.0,
+                          hole_points_UTM=[ 292110.784, 7676551.710 ])
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        number_of_volumes = fid.variables['volumes']
+        #print "number_of_volumes",len(number_of_volumes) 
+        assert num.allclose(16, len(number_of_volumes))
+        
+        fid.close()
+        self.delete_mux(files)
+        #print "sww_file", sww_file 
+        os.remove(sww_file)
+        
+    def test_urs_ungridded_holeII(self):
+
+        # Check that if using a hole that returns no triangles,
+        # urs_ungridded2sww removes the hole label.
+        
+        lat_long = [[-20.5, 114.5],
+                    [-20.6, 114.6],
+                    [-20.5, 115.5],
+                    [-20.6, 115.4],
+                    
+                    
+                    [-21.5, 114.5],
+                    [-21.4, 114.6],
+                    [-21.5, 115.5],
+                    [-21.4, 115.4]
+                    ]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        #Easting:  292110.784  Northing: 7676551.710 
+        #Latitude:   -21  0 ' 0.00000 ''  Longitude: 115  0 ' 0.00000 '' 
+
+        urs2sww(base_name, mean_stage=-240992.0,
+                          hole_points_UTM=[ 292110.784, 7676551.710 ])
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        fid = NetCDFFile(sww_file)
+        
+        volumes = fid.variables['volumes']
+        #print "number_of_volumes",len(volumes)
+
+        fid.close()
+        os.remove(sww_file)
+        
+        urs2sww(base_name, mean_stage=-240992.0)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        fid = NetCDFFile(sww_file)
+        
+        volumes_again = fid.variables['volumes']
+        #print "number_of_volumes",len(volumes_again) 
+        assert num.allclose(len(volumes_again),
+                            len(volumes))
+        fid.close()
+        os.remove(sww_file)
+        self.delete_mux(files) 
+        
+    def test_urs_ungridded2sww_mint_maxt (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
+                          mint=101, maxt=500)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = [0,100,200,300]
+       
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+        
+               #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+        # = n*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_x = n*(e+tide+n)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert num.allclose(answer_x, actual_x)  #Meters
+        
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+        # = e*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_y = -1*e*(e+tide+n)
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert num.allclose(answer_y, actual_y)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs_ungridded2sww_mint_maxtII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
+                          mint=0, maxt=100000)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, fid.variables['x'][:],
+                                                fid.variables['y'][:]))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = [0,100,200,300,400,500]
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        assert num.allclose(stage[0], x +tide)
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs_ungridded2sww_mint_maxtIII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        try:
+            urs2sww(base_name, mean_stage=tide,
+                          origin =(50,23432,4343),
+                          mint=301, maxt=399,
+                              verbose=self.verbose)
+        except: 
+            pass
+        else:
+            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
+
+        self.delete_mux(files)
+        
+    def test_urs_ungridded2sww_mint_maxt_bad (self):       
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        try:
+            urs2sww(base_name, mean_stage=tide,
+                          origin =(50,23432,4343),
+                          mint=301, maxt=301,
+                              verbose=self.verbose)
+        except: 
+            pass
+        else:
+            self.failUnless(0 ==1, 'Bad input did not throw exception error!')
+
+        self.delete_mux(files)
+
+        
+    def test_URS_points_needed_and_urs_ungridded2sww(self):
+        # This doesn't actually check anything
+        #  
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        time_step_count = 2
+        time_step = 400
+        tide = -200000
+        zone = 50
+
+        boundary_polygon = [[250000,7660000],[280000,7660000],
+                             [280000,7630000],[250000,7630000]]
+        geo=urs.calculate_boundary_points(boundary_polygon, zone,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount,
+                              verbose=self.verbose)
+        lat_long = geo.get_data_points(as_lat_long=True)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide,
+                          verbose=self.verbose)
+        self.delete_mux(files)
+        os.remove( base_name + '.sww')
+    
+    def cache_test_URS_points_needed_and_urs_ungridded2sww(self):
+        
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        time_step_count = 2
+        time_step = 400
+        tide = -200000
+        zone = 50
+
+        boundary_polygon = [[250000,7660000],[270000,7650000],
+                             [280000,7630000],[250000,7630000]]
+        geo=URS_points_needed(boundary_polygon, zone,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount, use_cache=True,
+                              verbose=True)
+        
+    def visual_test_URS_points_needed_and_urs_ungridded2sww(self):
+        
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        time_step_count = 2
+        time_step = 400
+        tide = -200000
+        zone = 50
+
+        boundary_polygon = [[250000,7660000],[270000,7650000],
+                             [280000,7630000],[250000,7630000]]
+        geo=URS_points_needed(boundary_polygon, zone,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount)
+        lat_long = geo.get_data_points(as_lat_long=True)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide)
+        self.delete_mux(files)
+        os.remove( base_name + '.sww')
+        # extend this so it interpolates onto the boundary.
+        # have it fail if there is NaN
+
+
+  
+    def test_urs_ungridded2swwII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 2
+        time_step = 400
+        tide = 9000000
+        geo_reference = Geo_reference(50, 3434543,34534543)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide, origin = geo_reference)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+        # = n*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_x = n*(e+tide+n)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert num.allclose(answer_x, actual_x)  #Meters
+        
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+        # = e*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_y = -1*e*(e+tide+n)
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert num.allclose(answer_y, actual_y)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+
+
+    def test_urs_ungridded2sww (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 2
+        time_step = 400
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs2sww(base_name, mean_stage=tide,
+                          verbose=self.verbose)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert num.allclose(ensure_numeric(times),
+                            ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+        # = n*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_x = n*(e+tide+n)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert num.allclose(answer_x, actual_x)  #Meters
+        
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+        # = e*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_y = -1*e*(e+tide+n)
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert num.allclose(answer_y, actual_y)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Dem2Pts,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py	(revision 7770)
@@ -1,9 +1,27 @@
+import numpy as num
+from Scientific.IO.NetCDF import NetCDFFile
+
+from anuga.file.urs import Read_urs
+
+from anuga.coordinate_transforms.redfearn import redfearn, \
+     convert_from_latlon_to_utm
+
+from anuga.geospatial_data.geospatial_data import ensure_absolute, \
+                                                    Geospatial_data
+
+from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
+                            
+from anuga.utilities.numerical_tools import ensure_numeric                            
+
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+
+from sww_file import Read_sww, Write_sww
+
+
 ################################################################################
 # CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
 ################################################################################
-
-WAVEHEIGHT_MUX_LABEL = '-z-mux'
-EAST_VELOCITY_LABEL =  '-e-mux'
-NORTH_VELOCITY_LABEL =  '-n-mux'
 
 ##
@@ -95,5 +113,5 @@
     mux = {}
     for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
+        mux[quantity] = Read_urs(file)
 
     # Could check that the depth is the same. (hashing)
@@ -143,5 +161,5 @@
     for i in range(a_mux.time_step_count):
         mux_times.append(a_mux.time_step * i)
-    (mux_times_start_i, mux_times_fin_i) = mux2sww_time(mux_times, mint, maxt)
+    (mux_times_start_i, mux_times_fin_i) = read_time_from_mux(mux_times, mint, maxt)
     times = mux_times[mux_times_start_i:mux_times_fin_i]
 
@@ -209,2 +227,24 @@
 
     # Do some conversions while writing the sww file
+
+
+def read_time_from_mux(mux_times, mint, maxt):
+    """
+        Read a list of mux times.
+        Return start and finish times which lie within the passed time period.
+    """
+
+    if mint == None:
+        mux_times_start_i = 0
+    else:
+        mux_times_start_i = num.searchsorted(mux_times, mint)
+
+    if maxt == None:
+        mux_times_fin_i = len(mux_times)
+    else:
+        maxt += 0.5 # so if you specify a time where there is
+                    # data that time will be included
+        mux_times_fin_i = num.searchsorted(mux_times, maxt)
+
+    return mux_times_start_i, mux_times_fin_i
+
Index: /trunk/anuga_core/source/anuga/file_conversion/urs2txt.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/urs2txt.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/file_conversion/urs2txt.py	(revision 7770)
@@ -1,3 +1,3 @@
-
+from anuga.file.urs import Read_urs
 
 ##
@@ -20,5 +20,5 @@
     mux = {}
     for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
+        mux[quantity] = Read_urs(file)
 
     # Could check that the depth is the same. (hashing)
Index: /trunk/anuga_core/source/anuga/pmesh/mesh_interface.py
===================================================================
--- /trunk/anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 7770)
@@ -26,4 +26,5 @@
                              interior_regions=None,
                              interior_holes=None,
+                             hole_tags=None,
                              poly_geo_reference=None,
                              mesh_geo_reference=None,
@@ -53,5 +54,7 @@
     throw an error
     
-    Interior_holes is a list of ploygons for each hole.
+    Interior_holes is a list of polygons for each hole.
+    hole_tags is an optional list of boundary tags for the holes, see
+                boundary_tags parameter.
 
     This function does not allow segments to share points - use underlying
@@ -92,4 +95,5 @@
               'interior_regions': interior_regions,
               'interior_holes': interior_holes,
+              'hole_tags': hole_tags,
               'poly_geo_reference': poly_geo_reference,
               'mesh_geo_reference': mesh_geo_reference,
@@ -127,4 +131,5 @@
                               interior_regions=None,
                               interior_holes=None,
+                              hole_tags=None,
                               poly_geo_reference=None,
                               mesh_geo_reference=None,
@@ -311,4 +316,5 @@
         for polygon in interior_holes:
             m.add_hole_from_polygon(polygon,
+                                    segment_tags=hole_tags,
                                     geo_reference=poly_geo_reference)
        
Index: /trunk/anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7770)
@@ -268,45 +268,4 @@
 
 
-##
-# @brief 
-# @param filename 
-# @param x 
-# @param y 
-# @param z 
-def write_obj(filename, x, y, z):
-    """Store x,y,z vectors into filename (obj format).
-
-       Vectors are assumed to have dimension (M,3) where
-       M corresponds to the number elements.
-       triangles are assumed to be disconnected
-
-       The three numbers in each vector correspond to three vertices,
-
-       e.g. the x coordinate of vertex 1 of element i is in x[i,1]
-    """
-
-    import os.path
-
-    root, ext = os.path.splitext(filename)
-    if ext == '.obj':
-        FN = filename
-    else:
-        FN = filename + '.obj'
-
-    outfile = open(FN, 'wb')
-    outfile.write("# Triangulation as an obj file\n")
-
-    M, N = x.shape
-    assert N == 3  #Assuming three vertices per element
-
-    for i in range(M):
-        for j in range(N):
-            outfile.write("v %f %f %f\n" % (x[i,j], y[i,j], z[i,j]))
-
-    for i in range(M):
-        base = i * N
-        outfile.write("f %d %d %d\n" % (base+1, base+2, base+3))
-
-    outfile.close()
 
 ##
@@ -949,495 +908,4 @@
         log.critical(msg)
 
-################################################################################
-# Functions to obtain diagnostics from sww files
-################################################################################
-
-
-
-##
-# @brief Get values for quantities interpolated to polyline midpoints from SWW.
-# @param filename Path to file to read.
-# @param quantity_names Quantity names to get.
-# @param polyline Representation of desired cross-section.
-# @param verbose True if this function is to be verbose.
-# @return (segments, i_func) where segments is a list of Triangle_intersection
-#         instances and i_func is an instance of Interpolation_function.
-# @note For 'polyline' assume absolute UTM coordinates.
-def get_interpolated_quantities_at_polyline_midpoints(filename,
-                                                      quantity_names=None,
-                                                      polyline=None,
-                                                      verbose=False):
-    """Get values for quantities interpolated to polyline midpoints from SWW
-
-    Input:
-        filename - Name of sww file
-        quantity_names - Names of quantities to load
-        polyline: Representation of desired cross section - it may contain
-                  multiple sections allowing for complex shapes. Assume
-                  absolute UTM coordinates.
-                  Format [[x0, y0], [x1, y1], ...]
-
-    Output:
-        segments: list of instances of class Triangle_intersection
-        interpolation_function: Instance of class Interpolation_function
-
-
-      This function is used by get_flow_through_cross_section and
-      get_energy_through_cross_section
-    """
-
-    from anuga.fit_interpolate.interpolate import Interpolation_function
-
-    # Get mesh and quantities from sww file
-    X = get_mesh_and_quantities_from_file(filename,
-                                          quantities=quantity_names,
-                                          verbose=verbose)
-    mesh, quantities, time = X
-
-    # Find all intersections and associated triangles.
-    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
-
-    # Get midpoints
-    interpolation_points = segment_midpoints(segments)
-
-    # Interpolate
-    if verbose:
-        log.critical('Interpolating - total number of interpolation points = %d'
-                     % len(interpolation_points))
-
-    I = Interpolation_function(time,
-                               quantities,
-                               quantity_names=quantity_names,
-                               vertex_coordinates=mesh.nodes,
-                               triangles=mesh.triangles,
-                               interpolation_points=interpolation_points,
-                               verbose=verbose)
-
-    return segments, I
-
-
-##
-# @brief Obtain flow (m^3/s) perpendicular to specified cross section.
-# @param filename Path to file to read.
-# @param polyline Representation of desired cross-section.
-# @param verbose Trie if this function is to be verbose.
-# @return (time, Q) where time and Q are lists of time and flow respectively.
-def get_flow_through_cross_section(filename, polyline, verbose=False):
-    """Obtain flow (m^3/s) perpendicular to specified cross section.
-
-    Inputs:
-        filename: Name of sww file
-        polyline: Representation of desired cross section - it may contain
-                  multiple sections allowing for complex shapes. Assume
-                  absolute UTM coordinates.
-                  Format [[x0, y0], [x1, y1], ...]
-
-    Output:
-        time: All stored times in sww file
-        Q: Hydrograph of total flow across given segments for all stored times.
-
-    The normal flow is computed for each triangle intersected by the polyline
-    and added up.  Multiple segments at different angles are specified the
-    normal flows may partially cancel each other.
-
-    The typical usage of this function would be to get flow through a channel,
-    and the polyline would then be a cross section perpendicular to the flow.
-    """
-
-    quantity_names =['elevation',
-                     'stage',
-                     'xmomentum',
-                     'ymomentum']
-
-    # Get values for quantities at each midpoint of poly line from sww file
-    X = get_interpolated_quantities_at_polyline_midpoints(filename,
-                                                          quantity_names=\
-                                                              quantity_names,
-                                                          polyline=polyline,
-                                                          verbose=verbose)
-    segments, interpolation_function = X
-
-    # Get vectors for time and interpolation_points
-    time = interpolation_function.time
-    interpolation_points = interpolation_function.interpolation_points
-
-    if verbose: log.critical('Computing hydrograph')
-
-    # Compute hydrograph
-    Q = []
-    for t in time:
-        total_flow = 0
-        for i in range(len(interpolation_points)):
-            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
-            normal = segments[i].normal
-
-            # Inner product of momentum vector with segment normal [m^2/s]
-            normal_momentum = uh*normal[0] + vh*normal[1]
-
-            # Flow across this segment [m^3/s]
-            segment_flow = normal_momentum * segments[i].length
-
-            # Accumulate
-            total_flow += segment_flow
-
-        # Store flow at this timestep
-        Q.append(total_flow)
-
-
-    return time, Q
-
-
-##
-# @brief Get average energy across a cross-section.
-# @param filename Path to file of interest.
-# @param polyline Representation of desired cross-section.
-# @param kind Select energy to compute: 'specific' or 'total'.
-# @param verbose True if this function is to be verbose.
-# @return (time, E) where time and E are lists of timestep and energy.
-def get_energy_through_cross_section(filename,
-                                     polyline,
-                                     kind='total',
-                                     verbose=False):
-    """Obtain average energy head [m] across specified cross section.
-
-    Inputs:
-        polyline: Representation of desired cross section - it may contain
-                  multiple sections allowing for complex shapes. Assume
-                  absolute UTM coordinates.
-                  Format [[x0, y0], [x1, y1], ...]
-        kind:     Select which energy to compute.
-                  Options are 'specific' and 'total' (default)
-
-    Output:
-        E: Average energy [m] across given segments for all stored times.
-
-    The average velocity is computed for each triangle intersected by
-    the polyline and averaged weighted by segment lengths.
-
-    The typical usage of this function would be to get average energy of
-    flow in a channel, and the polyline would then be a cross section
-    perpendicular to the flow.
-
-    #FIXME (Ole) - need name for this energy reflecting that its dimension
-    is [m].
-    """
-
-    from anuga.config import g, epsilon, velocity_protection as h0
-
-    quantity_names =['elevation',
-                     'stage',
-                     'xmomentum',
-                     'ymomentum']
-
-    # Get values for quantities at each midpoint of poly line from sww file
-    X = get_interpolated_quantities_at_polyline_midpoints(filename,
-                                                          quantity_names=\
-                                                              quantity_names,
-                                                          polyline=polyline,
-                                                          verbose=verbose)
-    segments, interpolation_function = X
-
-    # Get vectors for time and interpolation_points
-    time = interpolation_function.time
-    interpolation_points = interpolation_function.interpolation_points
-
-    if verbose: log.critical('Computing %s energy' % kind)
-
-    # Compute total length of polyline for use with weighted averages
-    total_line_length = 0.0
-    for segment in segments:
-        total_line_length += segment.length
-
-    # Compute energy
-    E = []
-    for t in time:
-        average_energy = 0.0
-        for i, p in enumerate(interpolation_points):
-            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
-
-            # Depth
-            h = depth = stage-elevation
-
-            # Average velocity across this segment
-            if h > epsilon:
-                # Use protection against degenerate velocities
-                u = uh / (h + h0/h)
-                v = vh / (h + h0/h)
-            else:
-                u = v = 0.0
-
-            speed_squared = u*u + v*v
-            kinetic_energy = 0.5 * speed_squared / g
-
-            if kind == 'specific':
-                segment_energy = depth + kinetic_energy
-            elif kind == 'total':
-                segment_energy = stage + kinetic_energy
-            else:
-                msg = 'Energy kind must be either "specific" or "total". '
-                msg += 'I got %s' % kind
-
-            # Add to weighted average
-            weigth = segments[i].length / total_line_length
-            average_energy += segment_energy * weigth
-
-        # Store energy at this timestep
-        E.append(average_energy)
-
-    return time, E
-
-
-##
-# @brief Return highest elevation where depth > 0.
-# @param filename Path to SWW file of interest.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param verbose True if this function is  to be verbose.
-def get_maximum_inundation_elevation(filename,
-                                     polygon=None,
-                                     time_interval=None,
-                                     verbose=False):
-    """Return highest elevation where depth > 0
-
-    Usage:
-    max_runup = get_maximum_inundation_elevation(filename,
-                                                 polygon=None,
-                                                 time_interval=None,
-                                                 verbose=False)
-
-    filename is a NetCDF sww file containing ANUGA model output.
-    Optional arguments polygon and time_interval restricts the maximum
-    runup calculation
-    to a points that lie within the specified polygon and time interval.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-
-    See general function get_maximum_inundation_data for details.
-    """
-
-    runup, _ = get_maximum_inundation_data(filename,
-                                           polygon=polygon,
-                                           time_interval=time_interval,
-                                           verbose=verbose)
-    return runup
-
-
-##
-# @brief Return location of highest elevation where h > 0
-# @param filename Path to SWW file to read.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param verbose True if this function is  to be verbose.
-def get_maximum_inundation_location(filename,
-                                    polygon=None,
-                                    time_interval=None,
-                                    verbose=False):
-    """Return location of highest elevation where h > 0
-
-    Usage:
-    max_runup_location = get_maximum_inundation_location(filename,
-                                                         polygon=None,
-                                                         time_interval=None,
-                                                         verbose=False)
-
-    filename is a NetCDF sww file containing ANUGA model output.
-    Optional arguments polygon and time_interval restricts the maximum
-    runup calculation
-    to a points that lie within the specified polygon and time interval.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-
-    See general function get_maximum_inundation_data for details.
-    """
-
-    _, max_loc = get_maximum_inundation_data(filename,
-                                             polygon=polygon,
-                                             time_interval=time_interval,
-                                             verbose=verbose)
-    return max_loc
-
-
-##
-# @brief Compute maximum run up height from SWW file.
-# @param filename Path to SWW file to read.
-# @param polygon If specified resrict to points inside this polygon.
-# @param time_interval If specified resrict to within the time specified.
-# @param use_centroid_values 
-# @param verbose True if this function is to be verbose.
-# @return (maximal_runup, maximal_runup_location)
-def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
-                                use_centroid_values=False,
-                                verbose=False):
-    """Compute maximum run up height from sww file.
-
-    Usage:
-    runup, location = get_maximum_inundation_data(filename,
-                                                  polygon=None,
-                                                  time_interval=None,
-                                                  verbose=False)
-
-    Algorithm is as in get_maximum_inundation_elevation from
-    shallow_water_domain except that this function works with the sww file and
-    computes the maximal runup height over multiple timesteps.
-
-    Optional arguments polygon and time_interval restricts the maximum runup
-    calculation to a points that lie within the specified polygon and time
-    interval.
-
-    Polygon is assumed to be in (absolute) UTM coordinates in the same zone
-    as domain.
-
-    If no inundation is found within polygon and time_interval the return value
-    is None signifying "No Runup" or "Everything is dry".
-    """
-
-    # We are using nodal values here as that is what is stored in sww files.
-
-    # Water depth below which it is considered to be 0 in the model
-    # FIXME (Ole): Allow this to be specified as a keyword argument as well
-
-    from anuga.geometry.polygon import inside_polygon
-    from anuga.config import minimum_allowed_height
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    dir, base = os.path.split(filename)
-
-    iterate_over = get_all_swwfiles(dir, base)
-
-    # Read sww file
-    if verbose: log.critical('Reading from %s' % filename)
-    # FIXME: Use general swwstats (when done)
-
-    maximal_runup = None
-    maximal_runup_location = None
-
-    for file, swwfile in enumerate (iterate_over):
-        # Read sww file
-        filename = join(dir, swwfile+'.sww')
-
-        if verbose: log.critical('Reading from %s' % filename)
-        # FIXME: Use general swwstats (when done)
-
-        fid = NetCDFFile(filename)
-
-        # 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()
-
-        # Get extent
-        volumes = fid.variables['volumes'][:]
-        x = fid.variables['x'][:] + xllcorner
-        y = fid.variables['y'][:] + yllcorner
-
-        # Get the relevant quantities (Convert from single precison)
-        elevation = num.array(fid.variables['elevation'][:], num.float)
-        stage = num.array(fid.variables['stage'][:], num.float)
-
-        # Here's where one could convert nodal information to centroid
-        # information but is probably something we need to write in C.
-        # Here's a Python thought which is NOT finished!!!
-        if use_centroid_values is True:
-            x = get_centroid_values(x, volumes)
-            y = get_centroid_values(y, volumes)
-            elevation = get_centroid_values(elevation, volumes)
-
-        # Spatial restriction
-        if polygon is not None:
-            msg = 'polygon must be a sequence of points.'
-            assert len(polygon[0]) == 2, msg
-            # FIXME (Ole): Make a generic polygon input check in polygon.py
-            # and call it here
-            points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
-                                                            y[:,num.newaxis]),
-                                                            axis=1))
-            point_indices = inside_polygon(points, polygon)
-
-            # Restrict quantities to polygon
-            elevation = num.take(elevation, point_indices, axis=0)
-            stage = num.take(stage, point_indices, axis=1)
-
-            # Get info for location of maximal runup
-            points_in_polygon = num.take(points, point_indices, axis=0)
-
-            x = points_in_polygon[:,0]
-            y = points_in_polygon[:,1]
-        else:
-            # Take all points
-            point_indices = num.arange(len(x))
-
-        # Temporal restriction
-        time = fid.variables['time'][:]
-        all_timeindices = num.arange(len(time))
-        if time_interval is not None:
-            msg = 'time_interval must be a sequence of length 2.'
-            assert len(time_interval) == 2, msg
-            msg = 'time_interval %s must not be decreasing.' % time_interval
-            assert time_interval[1] >= time_interval[0], msg
-            msg = 'Specified time interval [%.8f:%.8f] ' % tuple(time_interval)
-            msg += 'must does not match model time interval: [%.8f, %.8f]\n' \
-                   % (time[0], time[-1])
-            if time_interval[1] < time[0]: raise ValueError(msg)
-            if time_interval[0] > time[-1]: raise ValueError(msg)
-
-            # Take time indices corresponding to interval (& is bitwise AND)
-            timesteps = num.compress((time_interval[0] <= time) \
-                                     & (time <= time_interval[1]),
-                                     all_timeindices)
-
-            msg = 'time_interval %s did not include any model timesteps.' \
-                  % time_interval
-            assert not num.alltrue(timesteps == 0), msg
-        else:
-            # Take them all
-            timesteps = all_timeindices
-
-        fid.close()
-
-        # Compute maximal runup for each timestep
-        #maximal_runup = None
-        #maximal_runup_location = None
-        #maximal_runups = [None]
-        #maximal_runup_locations = [None]
-
-        for i in timesteps:
-            if use_centroid_values is True:
-                stage_i = get_centroid_values(stage[i,:], volumes)
-            else:
-                stage_i = stage[i,:]
-
-            depth = stage_i - elevation
-
-            # Get wet nodes i.e. nodes with depth>0 within given region
-            # and timesteps
-            wet_nodes = num.compress(depth > minimum_allowed_height,
-                                     num.arange(len(depth)))
-
-            if num.alltrue(wet_nodes == 0):
-                runup = None
-            else:
-                # Find maximum elevation among wet nodes
-                wet_elevation = num.take(elevation, wet_nodes, axis=0)
-                runup_index = num.argmax(wet_elevation)
-                runup = max(wet_elevation)
-                assert wet_elevation[runup_index] == runup       # Must be True
-
-            if runup > maximal_runup:
-                maximal_runup = runup      # works even if maximal_runup is None
-
-                # Record location
-                wet_x = num.take(x, wet_nodes, axis=0)
-                wet_y = num.take(y, wet_nodes, axis=0)
-                maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]
-
-    return maximal_runup, maximal_runup_location
 
 
Index: /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py	(revision 7770)
+++ /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py	(revision 7770)
@@ -0,0 +1,684 @@
+'''
+    Operations to extract information from an SWW file.
+'''
+
+##
+# @brief Get values for quantities interpolated to polyline midpoints from SWW.
+# @param filename Path to file to read.
+# @param quantity_names Quantity names to get.
+# @param polyline Representation of desired cross-section.
+# @param verbose True if this function is to be verbose.
+# @return (segments, i_func) where segments is a list of Triangle_intersection
+#         instances and i_func is an instance of Interpolation_function.
+# @note For 'polyline' assume absolute UTM coordinates.
+def get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                      quantity_names=None,
+                                                      polyline=None,
+                                                      verbose=False):
+    """Get values for quantities interpolated to polyline midpoints from SWW
+
+    Input:
+        filename - Name of sww file
+        quantity_names - Names of quantities to load
+        polyline: Representation of desired cross section - it may contain
+                  multiple sections allowing for complex shapes. Assume
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
+
+    Output:
+        segments: list of instances of class Triangle_intersection
+        interpolation_function: Instance of class Interpolation_function
+
+
+      This function is used by get_flow_through_cross_section and
+      get_energy_through_cross_section
+    """
+
+    from anuga.fit_interpolate.interpolate import Interpolation_function
+
+    # Get mesh and quantities from sww file
+    X = get_mesh_and_quantities_from_file(filename,
+                                          quantities=quantity_names,
+                                          verbose=verbose)
+    mesh, quantities, time = X
+
+    # Find all intersections and associated triangles.
+    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
+
+    # Get midpoints
+    interpolation_points = segment_midpoints(segments)
+
+    # Interpolate
+    if verbose:
+        log.critical('Interpolating - total number of interpolation points = %d'
+                     % len(interpolation_points))
+
+    I = Interpolation_function(time,
+                               quantities,
+                               quantity_names=quantity_names,
+                               vertex_coordinates=mesh.nodes,
+                               triangles=mesh.triangles,
+                               interpolation_points=interpolation_points,
+                               verbose=verbose)
+
+    return segments, I
+
+
+##
+# @brief Obtain flow (m^3/s) perpendicular to specified cross section.
+# @param filename Path to file to read.
+# @param polyline Representation of desired cross-section.
+# @param verbose Trie if this function is to be verbose.
+# @return (time, Q) where time and Q are lists of time and flow respectively.
+def get_flow_through_cross_section(filename, polyline, verbose=False):
+    """Obtain flow (m^3/s) perpendicular to specified cross section.
+
+    Inputs:
+        filename: Name of sww file
+        polyline: Representation of desired cross section - it may contain
+                  multiple sections allowing for complex shapes. Assume
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
+
+    Output:
+        time: All stored times in sww file
+        Q: Hydrograph of total flow across given segments for all stored times.
+
+    The normal flow is computed for each triangle intersected by the polyline
+    and added up.  Multiple segments at different angles are specified the
+    normal flows may partially cancel each other.
+
+    The typical usage of this function would be to get flow through a channel,
+    and the polyline would then be a cross section perpendicular to the flow.
+    """
+
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
+
+    # Get values for quantities at each midpoint of poly line from sww file
+    X = get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                          quantity_names=\
+                                                              quantity_names,
+                                                          polyline=polyline,
+                                                          verbose=verbose)
+    segments, interpolation_function = X
+
+    # Get vectors for time and interpolation_points
+    time = interpolation_function.time
+    interpolation_points = interpolation_function.interpolation_points
+
+    if verbose: log.critical('Computing hydrograph')
+
+    # Compute hydrograph
+    Q = []
+    for t in time:
+        total_flow = 0
+        for i in range(len(interpolation_points)):
+            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
+            normal = segments[i].normal
+
+            # Inner product of momentum vector with segment normal [m^2/s]
+            normal_momentum = uh*normal[0] + vh*normal[1]
+
+            # Flow across this segment [m^3/s]
+            segment_flow = normal_momentum * segments[i].length
+
+            # Accumulate
+            total_flow += segment_flow
+
+        # Store flow at this timestep
+        Q.append(total_flow)
+
+
+    return time, Q
+
+
+##
+# @brief Get average energy across a cross-section.
+# @param filename Path to file of interest.
+# @param polyline Representation of desired cross-section.
+# @param kind Select energy to compute: 'specific' or 'total'.
+# @param verbose True if this function is to be verbose.
+# @return (time, E) where time and E are lists of timestep and energy.
+def get_energy_through_cross_section(filename,
+                                     polyline,
+                                     kind='total',
+                                     verbose=False):
+    """Obtain average energy head [m] across specified cross section.
+
+    Inputs:
+        polyline: Representation of desired cross section - it may contain
+                  multiple sections allowing for complex shapes. Assume
+                  absolute UTM coordinates.
+                  Format [[x0, y0], [x1, y1], ...]
+        kind:     Select which energy to compute.
+                  Options are 'specific' and 'total' (default)
+
+    Output:
+        E: Average energy [m] across given segments for all stored times.
+
+    The average velocity is computed for each triangle intersected by
+    the polyline and averaged weighted by segment lengths.
+
+    The typical usage of this function would be to get average energy of
+    flow in a channel, and the polyline would then be a cross section
+    perpendicular to the flow.
+
+    #FIXME (Ole) - need name for this energy reflecting that its dimension
+    is [m].
+    """
+
+    from anuga.config import g, epsilon, velocity_protection as h0
+
+    quantity_names =['elevation',
+                     'stage',
+                     'xmomentum',
+                     'ymomentum']
+
+    # Get values for quantities at each midpoint of poly line from sww file
+    X = get_interpolated_quantities_at_polyline_midpoints(filename,
+                                                          quantity_names=\
+                                                              quantity_names,
+                                                          polyline=polyline,
+                                                          verbose=verbose)
+    segments, interpolation_function = X
+
+    # Get vectors for time and interpolation_points
+    time = interpolation_function.time
+    interpolation_points = interpolation_function.interpolation_points
+
+    if verbose: log.critical('Computing %s energy' % kind)
+
+    # Compute total length of polyline for use with weighted averages
+    total_line_length = 0.0
+    for segment in segments:
+        total_line_length += segment.length
+
+    # Compute energy
+    E = []
+    for t in time:
+        average_energy = 0.0
+        for i, p in enumerate(interpolation_points):
+            elevation, stage, uh, vh = interpolation_function(t, point_id=i)
+
+            # Depth
+            h = depth = stage-elevation
+
+            # Average velocity across this segment
+            if h > epsilon:
+                # Use protection against degenerate velocities
+                u = uh / (h + h0/h)
+                v = vh / (h + h0/h)
+            else:
+                u = v = 0.0
+
+            speed_squared = u*u + v*v
+            kinetic_energy = 0.5 * speed_squared / g
+
+            if kind == 'specific':
+                segment_energy = depth + kinetic_energy
+            elif kind == 'total':
+                segment_energy = stage + kinetic_energy
+            else:
+                msg = 'Energy kind must be either "specific" or "total". '
+                msg += 'I got %s' % kind
+
+            # Add to weighted average
+            weigth = segments[i].length / total_line_length
+            average_energy += segment_energy * weigth
+
+        # Store energy at this timestep
+        E.append(average_energy)
+
+    return time, E
+
+
+##
+# @brief Return highest elevation where depth > 0.
+# @param filename Path to SWW file of interest.
+# @param polygon If specified resrict to points inside this polygon.
+# @param time_interval If specified resrict to within the time specified.
+# @param verbose True if this function is  to be verbose.
+def get_maximum_inundation_elevation(filename,
+                                     polygon=None,
+                                     time_interval=None,
+                                     verbose=False):
+    """Return highest elevation where depth > 0
+
+    Usage:
+    max_runup = get_maximum_inundation_elevation(filename,
+                                                 polygon=None,
+                                                 time_interval=None,
+                                                 verbose=False)
+
+    filename is a NetCDF sww file containing ANUGA model output.
+    Optional arguments polygon and time_interval restricts the maximum
+    runup calculation
+    to a points that lie within the specified polygon and time interval.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
+
+    See general function get_maximum_inundation_data for details.
+    """
+
+    runup, _ = get_maximum_inundation_data(filename,
+                                           polygon=polygon,
+                                           time_interval=time_interval,
+                                           verbose=verbose)
+    return runup
+
+
+##
+# @brief Return location of highest elevation where h > 0
+# @param filename Path to SWW file to read.
+# @param polygon If specified resrict to points inside this polygon.
+# @param time_interval If specified resrict to within the time specified.
+# @param verbose True if this function is  to be verbose.
+def get_maximum_inundation_location(filename,
+                                    polygon=None,
+                                    time_interval=None,
+                                    verbose=False):
+    """Return location of highest elevation where h > 0
+
+    Usage:
+    max_runup_location = get_maximum_inundation_location(filename,
+                                                         polygon=None,
+                                                         time_interval=None,
+                                                         verbose=False)
+
+    filename is a NetCDF sww file containing ANUGA model output.
+    Optional arguments polygon and time_interval restricts the maximum
+    runup calculation
+    to a points that lie within the specified polygon and time interval.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
+
+    See general function get_maximum_inundation_data for details.
+    """
+
+    _, max_loc = get_maximum_inundation_data(filename,
+                                             polygon=polygon,
+                                             time_interval=time_interval,
+                                             verbose=verbose)
+    return max_loc
+
+
+##
+# @brief Compute maximum run up height from SWW file.
+# @param filename Path to SWW file to read.
+# @param polygon If specified resrict to points inside this polygon.
+# @param time_interval If specified resrict to within the time specified.
+# @param use_centroid_values 
+# @param verbose True if this function is to be verbose.
+# @return (maximal_runup, maximal_runup_location)
+def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
+                                use_centroid_values=False,
+                                verbose=False):
+    """Compute maximum run up height from sww file.
+
+    Usage:
+    runup, location = get_maximum_inundation_data(filename,
+                                                  polygon=None,
+                                                  time_interval=None,
+                                                  verbose=False)
+
+    Algorithm is as in get_maximum_inundation_elevation from
+    shallow_water_domain except that this function works with the sww file and
+    computes the maximal runup height over multiple timesteps.
+
+    Optional arguments polygon and time_interval restricts the maximum runup
+    calculation to a points that lie within the specified polygon and time
+    interval.
+
+    Polygon is assumed to be in (absolute) UTM coordinates in the same zone
+    as domain.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
+    """
+
+    # We are using nodal values here as that is what is stored in sww files.
+
+    # Water depth below which it is considered to be 0 in the model
+    # FIXME (Ole): Allow this to be specified as a keyword argument as well
+
+    from anuga.geometry.polygon import inside_polygon
+    from anuga.config import minimum_allowed_height
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    dir, base = os.path.split(filename)
+
+    iterate_over = get_all_swwfiles(dir, base)
+
+    # Read sww file
+    if verbose: log.critical('Reading from %s' % filename)
+    # FIXME: Use general swwstats (when done)
+
+    maximal_runup = None
+    maximal_runup_location = None
+
+    for file, swwfile in enumerate (iterate_over):
+        # Read sww file
+        filename = join(dir, swwfile+'.sww')
+
+        if verbose: log.critical('Reading from %s' % filename)
+        # FIXME: Use general swwstats (when done)
+
+        fid = NetCDFFile(filename)
+
+        # 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()
+
+        # Get extent
+        volumes = fid.variables['volumes'][:]
+        x = fid.variables['x'][:] + xllcorner
+        y = fid.variables['y'][:] + yllcorner
+
+        # Get the relevant quantities (Convert from single precison)
+        elevation = num.array(fid.variables['elevation'][:], num.float)
+        stage = num.array(fid.variables['stage'][:], num.float)
+
+        # Here's where one could convert nodal information to centroid
+        # information but is probably something we need to write in C.
+        # Here's a Python thought which is NOT finished!!!
+        if use_centroid_values is True:
+            x = get_centroid_values(x, volumes)
+            y = get_centroid_values(y, volumes)
+            elevation = get_centroid_values(elevation, volumes)
+
+        # Spatial restriction
+        if polygon is not None:
+            msg = 'polygon must be a sequence of points.'
+            assert len(polygon[0]) == 2, msg
+            # FIXME (Ole): Make a generic polygon input check in polygon.py
+            # and call it here
+            points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis],
+                                                            y[:,num.newaxis]),
+                                                            axis=1))
+            point_indices = inside_polygon(points, polygon)
+
+            # Restrict quantities to polygon
+            elevation = num.take(elevation, point_indices, axis=0)
+            stage = num.take(stage, point_indices, axis=1)
+
+            # Get info for location of maximal runup
+            points_in_polygon = num.take(points, point_indices, axis=0)
+
+            x = points_in_polygon[:,0]
+            y = points_in_polygon[:,1]
+        else:
+            # Take all points
+            point_indices = num.arange(len(x))
+
+        # Temporal restriction
+        time = fid.variables['time'][:]
+        all_timeindices = num.arange(len(time))
+        if time_interval is not None:
+            msg = 'time_interval must be a sequence of length 2.'
+            assert len(time_interval) == 2, msg
+            msg = 'time_interval %s must not be decreasing.' % time_interval
+            assert time_interval[1] >= time_interval[0], msg
+            msg = 'Specified time interval [%.8f:%.8f] ' % tuple(time_interval)
+            msg += 'must does not match model time interval: [%.8f, %.8f]\n' \
+                   % (time[0], time[-1])
+            if time_interval[1] < time[0]: raise ValueError(msg)
+            if time_interval[0] > time[-1]: raise ValueError(msg)
+
+            # Take time indices corresponding to interval (& is bitwise AND)
+            timesteps = num.compress((time_interval[0] <= time) \
+                                     & (time <= time_interval[1]),
+                                     all_timeindices)
+
+            msg = 'time_interval %s did not include any model timesteps.' \
+                  % time_interval
+            assert not num.alltrue(timesteps == 0), msg
+        else:
+            # Take them all
+            timesteps = all_timeindices
+
+        fid.close()
+
+        # Compute maximal runup for each timestep
+        #maximal_runup = None
+        #maximal_runup_location = None
+        #maximal_runups = [None]
+        #maximal_runup_locations = [None]
+
+        for i in timesteps:
+            if use_centroid_values is True:
+                stage_i = get_centroid_values(stage[i,:], volumes)
+            else:
+                stage_i = stage[i,:]
+
+            depth = stage_i - elevation
+
+            # Get wet nodes i.e. nodes with depth>0 within given region
+            # and timesteps
+            wet_nodes = num.compress(depth > minimum_allowed_height,
+                                     num.arange(len(depth)))
+
+            if num.alltrue(wet_nodes == 0):
+                runup = None
+            else:
+                # Find maximum elevation among wet nodes
+                wet_elevation = num.take(elevation, wet_nodes, axis=0)
+                runup_index = num.argmax(wet_elevation)
+                runup = max(wet_elevation)
+                assert wet_elevation[runup_index] == runup       # Must be True
+
+            if runup > maximal_runup:
+                maximal_runup = runup      # works even if maximal_runup is None
+
+                # Record location
+                wet_x = num.take(x, wet_nodes, axis=0)
+                wet_y = num.take(y, wet_nodes, axis=0)
+                maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]
+
+    return maximal_runup, maximal_runup_location
+
+
+
+    def test_get_maximum_inundation_from_sww(self):
+        """test_get_maximum_inundation_from_sww(self)
+
+        Test of get_maximum_inundation_elevation()
+        and get_maximum_inundation_location() from data_manager.py
+
+        This is based on test_get_maximum_inundation_3(self) but works with the
+        stored results instead of with the internal data structure.
+
+        This test uses the underlying get_maximum_inundation_data for tests
+        """
+
+        from anuga.abstract_2d_finite_volumes.mesh_factory \
+                import rectangular_cross
+        from data_manager import get_maximum_inundation_elevation
+        from data_manager import get_maximum_inundation_location
+        from data_manager import get_maximum_inundation_data
+
+        initial_runup_height = -0.4
+        final_runup_height = -0.3
+
+        #--------------------------------------------------------------
+        # Setup computational domain
+        #--------------------------------------------------------------
+        N = 10
+        points, vertices, boundary = rectangular_cross(N, N)
+        domain = Domain(points, vertices, boundary)
+        domain.set_name('runup_test')
+        domain.set_maximum_allowed_speed(1.0)
+
+        # FIXME: This works better with old limiters so far
+        domain.tight_slope_limiters = 0
+
+        #--------------------------------------------------------------
+        # Setup initial conditions
+        #--------------------------------------------------------------
+        def topography(x, y):
+            return -x/2                             # linear bed slope
+
+        # Use function for elevation
+        domain.set_quantity('elevation', topography)
+        domain.set_quantity('friction', 0.)                # Zero friction
+        # Constant negative initial stage
+        domain.set_quantity('stage', initial_runup_height)
+
+        #--------------------------------------------------------------
+        # Setup boundary conditions
+        #--------------------------------------------------------------
+        Br = Reflective_boundary(domain)                       # Reflective wall
+        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
+
+        # All reflective to begin with (still water)
+        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+        #--------------------------------------------------------------
+        # Test initial inundation height
+        #--------------------------------------------------------------
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+                get_values(location='centroids', indices=indices)
+        assert num.alltrue(z < initial_runup_height)
+
+        q_ref = domain.get_maximum_inundation_elevation()
+        # First order accuracy
+        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
+
+        #--------------------------------------------------------------
+        # Let triangles adjust
+        #--------------------------------------------------------------
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
+            pass
+
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+        q_ref = domain.get_maximum_inundation_elevation()
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        msg = 'We got %f, should have been %f' % (q, q_ref)
+        assert num.allclose(q, q_ref, rtol=1.0/N), msg
+
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
+        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
+
+        # Test error condition if time interval is out
+        try:
+            q = get_maximum_inundation_elevation('runup_test.sww',
+                                                 time_interval=[2.0, 3.0])
+        except ValueError:
+            pass
+        else:
+            msg = 'should have caught wrong time interval'
+            raise Exception, msg
+
+        # Check correct time interval
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             time_interval=[0.0, 3.0])
+        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
+        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
+        assert num.allclose(-loc[0]/2, q)    # From topography formula
+
+        #--------------------------------------------------------------
+        # Update boundary to allow inflow
+        #--------------------------------------------------------------
+        domain.set_boundary({'right': Bd})
+
+        #--------------------------------------------------------------
+        # Evolve system through time
+        #--------------------------------------------------------------
+        q_max = None
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
+                               skip_initial_step = True):
+            q = domain.get_maximum_inundation_elevation()
+            if q > q_max:
+                q_max = q
+
+        #--------------------------------------------------------------
+        # Test inundation height again
+        #--------------------------------------------------------------
+        indices = domain.get_wet_elements()
+        z = domain.get_quantity('elevation').\
+                get_values(location='centroids', indices=indices)
+
+        assert num.alltrue(z < final_runup_height)
+
+        q = domain.get_maximum_inundation_elevation()
+        # First order accuracy
+        assert num.allclose(q, final_runup_height, rtol=1.0/N)
+
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             time_interval=[3.0, 3.0])
+        msg = 'We got %f, should have been %f' % (q, final_runup_height)
+        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
+        assert num.allclose(-loc[0]/2, q)    # From topography formula
+
+        q = get_maximum_inundation_elevation('runup_test.sww')
+        loc = get_maximum_inundation_location('runup_test.sww')
+        msg = 'We got %f, should have been %f' % (q, q_max)
+        assert num.allclose(q, q_max, rtol=1.0/N), msg
+        assert num.allclose(-loc[0]/2, q)    # From topography formula
+
+        q = get_maximum_inundation_elevation('runup_test.sww',
+                                             time_interval=[0, 3])
+        msg = 'We got %f, should have been %f' % (q, q_max)
+        assert num.allclose(q, q_max, rtol=1.0/N), msg
+
+        # Check polygon mode
+        # Runup region
+        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
+        q = get_maximum_inundation_elevation('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %f, should have been %f' % (q, q_max)
+        assert num.allclose(q, q_max, rtol=1.0/N), msg
+
+        # Offshore region
+        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %f, should have been %f' % (q, -0.475)
+        assert num.allclose(q, -0.475, rtol=1.0/N), msg
+        assert is_inside_polygon(loc, polygon)
+        assert num.allclose(-loc[0]/2, q)    # From topography formula
+
+        # Dry region
+        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %s, should have been None' % (q)
+        assert q is None, msg
+        msg = 'We got %s, should have been None' % (loc)
+        assert loc is None, msg
+
+        # Check what happens if no time point is within interval
+        try:
+            q = get_maximum_inundation_elevation('runup_test.sww',
+                                                 time_interval=[2.75, 2.75])
+        except AssertionError:
+            pass
+        else:
+            msg = 'Time interval should have raised an exception'
+            raise Exception, msg
+
+        # Cleanup
+        try:
+            os.remove(domain.get_name() + '.sww')
+        except:
+            pass
+            #FIXME(Ole): Windows won't allow removal of this
+
Index: /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7770)
@@ -1351,26 +1351,4 @@
         assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
 
-
-
-    def test_weed(self):
-        from data_manager import weed
-
-        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
-        volumes1 = [[0,1,2],[3,4,5]]
-        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
-        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
-
-        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
-
-        assert len(points2)==len(coordinates2)
-        for i in range(len(coordinates2)):
-            coordinate = tuple(coordinates2[i])
-            assert points2.has_key(coordinate)
-            points2[coordinate]=i
-
-        for triangle in volumes1:
-            for coordinate in triangle:
-                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
-                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
 
 
@@ -2488,552 +2466,4 @@
         
 
-    def test_get_maximum_inundation(self):
-        """Test that sww information can be converted correctly to maximum
-        runup elevation and location (without and with georeferencing)
-
-        This test creates a slope and a runup which is maximal (~11m) at around 10s
-        and levels out to the boundary condition (1m) at about 30s.
-        """
-
-        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(20, 5, 100, 50)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('runuptest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        # FIXME (Ole): Backwards compatibility
-        # Look at sww file and see what happens when
-        # domain.tight_slope_limiters = 1
-        domain.tight_slope_limiters = 0
-        domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)        
-        
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([1.0,0,0])
-
-
-        #---------- First run without geo referencing
-        
-        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
-        domain.set_quantity('stage', -6)
-        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = 50):
-            pass
-
-
-        # Check maximal runup
-        runup = get_maximum_inundation_elevation(swwfile)
-        location = get_maximum_inundation_location(swwfile)
-        #print 'Runup, location', runup, location
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15) or num.allclose(location[0], 10)
-
-        # Check final runup
-        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
-        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
-        # print 'Runup, location:',runup, location        
-        assert num.allclose(runup, 1)
-        assert num.allclose(location[0], 65)
-
-        # Check runup restricted to a polygon
-        p = [[50,1], [99,1], [99,49], [50,49]]
-        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
-        location = get_maximum_inundation_location(swwfile, polygon=p)
-        #print runup, location        
-        assert num.allclose(runup, 4)
-        assert num.allclose(location[0], 50)                
-
-        # Check that mimimum_storable_height works
-        fid = NetCDFFile(swwfile, netcdf_mode_r) # Open existing file
-        
-        stage = fid.variables['stage'][:]
-        z = fid.variables['elevation'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]        
-
-        
-        
-        for i in range(stage.shape[0]):
-            h = stage[i]-z # depth vector at time step i
-            
-            # Check every node location
-            for j in range(stage.shape[1]):
-                # Depth being either exactly zero implies
-                # momentum being zero.
-                # Or else depth must be greater than or equal to
-                # the minimal storable height
-                if h[j] == 0.0:
-                    assert xmomentum[i,j] == 0.0
-                    assert ymomentum[i,j] == 0.0                
-                else:
-                    assert h[j] >= domain.minimum_storable_height
-        
-        fid.close()
-
-        # Cleanup
-        os.remove(swwfile)
-        
-
-
-        #------------- Now the same with georeferencing
-
-        domain.time=0.0
-        E = 308500
-        N = 6189000
-        #E = N = 0
-        domain.geo_reference = Geo_reference(56, E, N)
-
-        domain.set_quantity('elevation', lambda x,y: -0.2*x + 14) # Slope
-        domain.set_quantity('stage', -6)
-        domain.set_boundary( {'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = 50):
-            pass
-
-        # Check maximal runup
-        runup = get_maximum_inundation_elevation(swwfile)
-        location = get_maximum_inundation_location(swwfile)
-        assert num.allclose(runup, 11) or num.allclose(runup, 12) # old limiters
-        assert num.allclose(location[0], 15+E) or num.allclose(location[0], 10+E)
-
-        # Check final runup
-        runup = get_maximum_inundation_elevation(swwfile, time_interval=[45,50])
-        location = get_maximum_inundation_location(swwfile, time_interval=[45,50])
-        assert num.allclose(runup, 1)
-        assert num.allclose(location[0], 65+E)
-
-        # Check runup restricted to a polygon
-        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.int) + num.array([E, N], num.int)      #array default#
-
-        runup = get_maximum_inundation_elevation(swwfile, polygon=p)
-        location = get_maximum_inundation_location(swwfile, polygon=p)
-        assert num.allclose(runup, 4)
-        assert num.allclose(location[0], 50+E)                
-
-
-        # Cleanup
-        os.remove(swwfile)
-
-
-    def test_get_mesh_and_quantities_from_sww_file(self):
-        """test_get_mesh_and_quantities_from_sww_file(self):
-        """     
-        
-        # Generate a test sww file with non trivial georeference
-        
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (100m x 5m)
-        width = 5
-        length = 50
-        t_end = 10
-        points, vertices, boundary = rectangular(length, width, 50, 5)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-        domain.set_datadir('.')
-
-        Br = Reflective_boundary(domain)    # Side walls
-        Bd = Dirichlet_boundary([1, 0, 0])  # inflow
-
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        
-        # Read it
-
-        # Get mesh and quantities from sww file
-        X = get_mesh_and_quantities_from_file(swwfile,
-                                              quantities=['elevation',
-                                                          'stage',
-                                                          'xmomentum',
-                                                          'ymomentum'], 
-                                              verbose=False)
-        mesh, quantities, time = X
-        
-
-        # Check that mesh has been recovered
-        assert num.alltrue(mesh.triangles == domain.get_triangles())
-        assert num.allclose(mesh.nodes, domain.get_nodes())
-
-        # Check that time has been recovered
-        assert num.allclose(time, range(t_end+1))
-
-        # Check that quantities have been recovered
-        # (sww files use single precision)
-        z=domain.get_quantity('elevation').get_values(location='unique vertices')
-        assert num.allclose(quantities['elevation'], z)
-
-        for q in ['stage', 'xmomentum', 'ymomentum']:
-            # Get quantity at last timestep
-            q_ref=domain.get_quantity(q).get_values(location='unique vertices')
-
-            #print q,quantities[q]
-            q_sww=quantities[q][-1,:]
-
-            msg = 'Quantity %s failed to be recovered' %q
-            assert num.allclose(q_ref, q_sww, atol=1.0e-6), msg
-            
-        # Cleanup
-        os.remove(swwfile)
-        
-
-    def test_get_flow_through_cross_section(self):
-        """test_get_flow_through_cross_section(self):
-
-        Test that the total flow through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected flow.
-
-        The specifics are
-        u = 2 m/s
-        h = 1 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 6 m^3/s
-
-        #---------- First run without geo referencing        
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 3
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        h = 1.0
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([h, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-
-        
-        domain.set_quantity('elevation', 0.0)
-        domain.set_quantity('stage', h)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-        for t in range(t_end+1):
-            for i in range(3):
-                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
-            
-
-        # Check flows through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            assert num.allclose(Q, uh*width)
-
-
-       
-        # Try the same with partial lines
-        x = length/2.
-        for i in range(5):
-            start_point = [length/2., i*width/5.]
-            #print start_point
-                            
-            cross_section = [start_point, [length/2., width]]
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            #print i, Q, (width-start_point[1])
-            assert num.allclose(Q, uh*(width-start_point[1]))
-
-
-        # Verify no flow when line is parallel to flow
-        cross_section = [[length/2.-10, width/2.], [length/2.+10, width/2.]]
-        time, Q = get_flow_through_cross_section(swwfile,
-                                                 cross_section,
-                                                 verbose=False)
-
-        #print i, Q
-        assert num.allclose(Q, 0, atol=1.0e-5)        
-
-
-        # Try with lines on an angle (all flow still runs through here)
-        cross_section = [[length/2., 0], [length/2.+width, width]]
-        time, Q = get_flow_through_cross_section(swwfile,
-                                                 cross_section,
-                                                 verbose=False)
-
-        assert num.allclose(Q, uh*width)        
-        
-
-
-                                      
-    def test_get_flow_through_cross_section_with_geo(self):
-        """test_get_flow_through_cross_section(self):
-
-        Test that the total flow through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected flow.
-
-        The specifics are
-        u = 2 m/s
-        h = 2 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 12 m^3/s
-
-
-        This run tries it with georeferencing and with elevation = -1
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 1
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        e = -1.0
-        w = 1.0
-        h = w-e
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-
-
-        
-        domain.set_quantity('elevation', e)
-        domain.set_quantity('stage', w)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        I = domain.geo_reference.get_absolute(I)
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-
-        for t in range(t_end+1):
-            for i in range(3):
-                #print i, t, f(t, i)            
-                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
-            
-
-        # Check flows through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-
-            cross_section = domain.geo_reference.get_absolute(cross_section)            
-            time, Q = get_flow_through_cross_section(swwfile,
-                                                     cross_section,
-                                                     verbose=False)
-
-            assert num.allclose(Q, uh*width)
-
-
-            
-    def test_get_energy_through_cross_section(self):
-        """test_get_energy_through_cross_section(self):
-
-        Test that the specific and total energy through a cross section can be
-        correctly obtained from an sww file.
-        
-        This test creates a flat bed with a known flow through it and tests
-        that the function correctly returns the expected energies.
-
-        The specifics are
-        u = 2 m/s
-        h = 1 m
-        w = 3 m (width of channel)
-
-        q = u*h*w = 6 m^3/s
-        Es = h + 0.5*v*v/g  # Specific energy head [m]
-        Et = w + 0.5*v*v/g  # Total energy head [m]        
-
-
-        This test uses georeferencing
-        
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        # Setup
-        from mesh_factory import rectangular
-
-        # Create basic mesh (20m x 3m)
-        width = 3
-        length = 20
-        t_end = 1
-        points, vertices, boundary = rectangular(length, width,
-                                                 length, width)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary,
-                        geo_reference = Geo_reference(56,308500,6189000))
-
-        domain.default_order = 2
-        domain.set_minimum_storable_height(0.01)
-
-        domain.set_name('flowtest')
-        swwfile = domain.get_name() + '.sww'
-
-        domain.set_datadir('.')
-        domain.format = 'sww'
-        domain.smooth = True
-
-        e = -1.0
-        w = 1.0
-        h = w-e
-        u = 2.0
-        uh = u*h
-
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
-
-        
-        domain.set_quantity('elevation', e)
-        domain.set_quantity('stage', w)
-        domain.set_quantity('xmomentum', uh)
-        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
-
-        for t in domain.evolve(yieldstep=1, finaltime = t_end):
-            pass
-
-        # Check that momentum is as it should be in the interior
-
-        I = [[0, width/2.],
-             [length/2., width/2.],
-             [length, width/2.]]
-        
-        I = domain.geo_reference.get_absolute(I)
-        f = file_function(swwfile,
-                          quantities=['stage', 'xmomentum', 'ymomentum'],
-                          interpolation_points=I,
-                          verbose=False)
-
-        for t in range(t_end+1):
-            for i in range(3):
-                #print i, t, f(t, i)
-                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
-            
-
-        # Check energies through the middle
-        for i in range(5):
-            x = length/2. + i*0.23674563 # Arbitrary
-            cross_section = [[x, 0], [x, width]]
-
-            cross_section = domain.geo_reference.get_absolute(cross_section)            
-            
-            time, Es = get_energy_through_cross_section(swwfile,
-                                                       cross_section,
-                                                       kind='specific',
-                                                       verbose=False)
-            assert num.allclose(Es, h + 0.5*u*u/g)
-            
-            time, Et = get_energy_through_cross_section(swwfile,
-                                                        cross_section,
-                                                        kind='total',
-                                                        verbose=False)
-            assert num.allclose(Et, w + 0.5*u*u/g)            
-
-            
         
     def test_get_all_swwfiles(self):
Index: /trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7769)
+++ /trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7770)
@@ -29,5 +29,5 @@
 from anuga.shallow_water.forcing import Rainfall, Wind_stress
 from anuga.shallow_water.forcing import Inflow, Cross_section
-from anuga.shallow_water.data_manager import get_flow_through_cross_section
+from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
 
 # boundary functions
@@ -1257,196 +1257,4 @@
         assert num.alltrue(wet_elevations < final_runup_height)
         assert num.allclose(wet_elevations, z)
-
-    def test_get_maximum_inundation_from_sww(self):
-        """test_get_maximum_inundation_from_sww(self)
-
-        Test of get_maximum_inundation_elevation()
-        and get_maximum_inundation_location() from data_manager.py
-
-        This is based on test_get_maximum_inundation_3(self) but works with the
-        stored results instead of with the internal data structure.
-
-        This test uses the underlying get_maximum_inundation_data for tests
-        """
-
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from data_manager import get_maximum_inundation_elevation
-        from data_manager import get_maximum_inundation_location
-        from data_manager import get_maximum_inundation_data
-
-        initial_runup_height = -0.4
-        final_runup_height = -0.3
-
-        #--------------------------------------------------------------
-        # Setup computational domain
-        #--------------------------------------------------------------
-        N = 10
-        points, vertices, boundary = rectangular_cross(N, N)
-        domain = Domain(points, vertices, boundary)
-        domain.set_name('runup_test')
-        domain.set_maximum_allowed_speed(1.0)
-
-        # FIXME: This works better with old limiters so far
-        domain.tight_slope_limiters = 0
-
-        #--------------------------------------------------------------
-        # Setup initial conditions
-        #--------------------------------------------------------------
-        def topography(x, y):
-            return -x/2                             # linear bed slope
-
-        # Use function for elevation
-        domain.set_quantity('elevation', topography)
-        domain.set_quantity('friction', 0.)                # Zero friction
-        # Constant negative initial stage
-        domain.set_quantity('stage', initial_runup_height)
-
-        #--------------------------------------------------------------
-        # Setup boundary conditions
-        #--------------------------------------------------------------
-        Br = Reflective_boundary(domain)                       # Reflective wall
-        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
-
-        # All reflective to begin with (still water)
-        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
-
-        #--------------------------------------------------------------
-        # Test initial inundation height
-        #--------------------------------------------------------------
-        indices = domain.get_wet_elements()
-        z = domain.get_quantity('elevation').\
-                get_values(location='centroids', indices=indices)
-        assert num.alltrue(z < initial_runup_height)
-
-        q_ref = domain.get_maximum_inundation_elevation()
-        # First order accuracy
-        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
-
-        #--------------------------------------------------------------
-        # Let triangles adjust
-        #--------------------------------------------------------------
-        for t in domain.evolve(yieldstep = 0.1, finaltime = 1.0):
-            pass
-
-        #--------------------------------------------------------------
-        # Test inundation height again
-        #--------------------------------------------------------------
-        q_ref = domain.get_maximum_inundation_elevation()
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, q_ref)
-        assert num.allclose(q, q_ref, rtol=1.0/N), msg
-
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
-        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
-
-        # Test error condition if time interval is out
-        try:
-            q = get_maximum_inundation_elevation('runup_test.sww',
-                                                 time_interval=[2.0, 3.0])
-        except ValueError:
-            pass
-        else:
-            msg = 'should have caught wrong time interval'
-            raise Exception, msg
-
-        # Check correct time interval
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             time_interval=[0.0, 3.0])
-        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
-        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        #--------------------------------------------------------------
-        # Update boundary to allow inflow
-        #--------------------------------------------------------------
-        domain.set_boundary({'right': Bd})
-
-        #--------------------------------------------------------------
-        # Evolve system through time
-        #--------------------------------------------------------------
-        q_max = None
-        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
-                               skip_initial_step = True):
-            q = domain.get_maximum_inundation_elevation()
-            if q > q_max:
-                q_max = q
-
-        #--------------------------------------------------------------
-        # Test inundation height again
-        #--------------------------------------------------------------
-        indices = domain.get_wet_elements()
-        z = domain.get_quantity('elevation').\
-                get_values(location='centroids', indices=indices)
-
-        assert num.alltrue(z < final_runup_height)
-
-        q = domain.get_maximum_inundation_elevation()
-        # First order accuracy
-        assert num.allclose(q, final_runup_height, rtol=1.0/N)
-
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             time_interval=[3.0, 3.0])
-        msg = 'We got %f, should have been %f' % (q, final_runup_height)
-        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        q = get_maximum_inundation_elevation('runup_test.sww')
-        loc = get_maximum_inundation_location('runup_test.sww')
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        q = get_maximum_inundation_elevation('runup_test.sww',
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-
-        # Check polygon mode
-        # Runup region
-        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
-        q = get_maximum_inundation_elevation('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, q_max)
-        assert num.allclose(q, q_max, rtol=1.0/N), msg
-
-        # Offshore region
-        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %f, should have been %f' % (q, -0.475)
-        assert num.allclose(q, -0.475, rtol=1.0/N), msg
-        assert is_inside_polygon(loc, polygon)
-        assert num.allclose(-loc[0]/2, q)    # From topography formula
-
-        # Dry region
-        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
-        q, loc = get_maximum_inundation_data('runup_test.sww',
-                                             polygon = polygon,
-                                             time_interval=[0, 3])
-        msg = 'We got %s, should have been None' % (q)
-        assert q is None, msg
-        msg = 'We got %s, should have been None' % (loc)
-        assert loc is None, msg
-
-        # Check what happens if no time point is within interval
-        try:
-            q = get_maximum_inundation_elevation('runup_test.sww',
-                                                 time_interval=[2.75, 2.75])
-        except AssertionError:
-            pass
-        else:
-            msg = 'Time interval should have raised an exception'
-            raise Exception, msg
-
-        # Cleanup
-        try:
-            os.remove(domain.get_name() + '.sww')
-        except:
-            pass
-            #FIXME(Ole): Windows won't allow removal of this
 
     def test_get_flow_through_cross_section_with_geo(self):
@@ -6085,5 +5893,5 @@
         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
         from anuga.shallow_water.shallow_water_domain import Rainfall
-        from anuga.shallow_water.data_manager import get_flow_through_cross_section
+        from anuga.shallow_water.sww_interrogate import get_flow_through_cross_section
 
         #----------------------------------------------------------------------
