Index: /trunk/anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7768)
+++ /trunk/anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7769)
@@ -455,326 +455,4 @@
 
 
-
-##
-# @brief Get the extents of a NetCDF data file.
-# @param file_name The path to the SWW file.
-# @return A list of x, y, z and stage limits (min, max).
-def extent_sww(file_name):
-    """Read in an sww file.
-
-    Input:
-    file_name - the sww file
-
-    Output:
-    A list: [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
-    """
-
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    #Get NetCDF
-    fid = NetCDFFile(file_name, netcdf_mode_r)
-
-    # Get the variables
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    stage = fid.variables['stage'][:]
-
-    fid.close()
-
-    return [min(x), max(x), min(y), max(y), num.min(stage), num.max(stage)]
-
-
-##
-# @brief 
-# @param filename
-# @param boundary
-# @param t
-# @param fail_if_NaN
-# @param NaN_filler
-# @param verbose
-# @param very_verbose
-# @return 
-def load_sww_as_domain(filename, boundary=None, t=None,
-               fail_if_NaN=True, NaN_filler=0,
-               verbose=False, very_verbose=False):
-    """
-    Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
-
-    Boundary is not recommended if domain.smooth is not selected, as it
-    uses unique coordinates, but not unique boundaries. This means that
-    the boundary file will not be compatable with the coordinates, and will
-    give a different final boundary, or crash.
-    """
-    
-    from Scientific.IO.NetCDF import NetCDFFile
-    from shallow_water import Domain
-
-    # initialise NaN.
-    NaN = 9.969209968386869e+036
-
-    if verbose: log.critical('Reading from %s' % filename)
-
-    fid = NetCDFFile(filename, netcdf_mode_r)    # Open existing file for read
-    time = fid.variables['time']       # Timesteps
-    if t is None:
-        t = time[-1]
-    time_interp = get_time_interp(time,t)
-
-    # Get the variables as numeric arrays
-    x = fid.variables['x'][:]                   # x-coordinates of vertices
-    y = fid.variables['y'][:]                   # y-coordinates of vertices
-    elevation = fid.variables['elevation']      # Elevation
-    stage = fid.variables['stage']              # Water level
-    xmomentum = fid.variables['xmomentum']      # Momentum in the x-direction
-    ymomentum = fid.variables['ymomentum']      # Momentum in the y-direction
-
-    starttime = fid.starttime[0]
-    volumes = fid.variables['volumes'][:]       # Connectivity
-    coordinates = num.transpose(num.asarray([x.tolist(), y.tolist()]))
-    # FIXME (Ole): Something like this might be better:
-    #                 concatenate((x, y), axis=1)
-    # or              concatenate((x[:,num.newaxis], x[:,num.newaxis]), axis=1)
-
-    conserved_quantities = []
-    interpolated_quantities = {}
-    other_quantities = []
-
-    # get geo_reference
-    try:                             # sww files don't have to have a geo_ref
-        geo_reference = Geo_reference(NetCDFObject=fid)
-    except: # AttributeError, e:
-        geo_reference = None
-
-    if verbose: log.critical('    getting quantities')
-
-    for quantity in fid.variables.keys():
-        dimensions = fid.variables[quantity].dimensions
-        if 'number_of_timesteps' in dimensions:
-            conserved_quantities.append(quantity)
-            interpolated_quantities[quantity] = \
-                  interpolated_quantity(fid.variables[quantity][:], time_interp)
-        else:
-            other_quantities.append(quantity)
-
-    other_quantities.remove('x')
-    other_quantities.remove('y')
-    #other_quantities.remove('z')
-    other_quantities.remove('volumes')
-    try:
-        other_quantities.remove('stage_range')
-        other_quantities.remove('xmomentum_range')
-        other_quantities.remove('ymomentum_range')
-        other_quantities.remove('elevation_range')
-    except:
-        pass
-
-    conserved_quantities.remove('time')
-
-    if verbose: log.critical('    building domain')
-
-    #    From domain.Domain:
-    #    domain = Domain(coordinates, volumes,\
-    #                    conserved_quantities = conserved_quantities,\
-    #                    other_quantities = other_quantities,zone=zone,\
-    #                    xllcorner=xllcorner, yllcorner=yllcorner)
-
-    # From shallow_water.Domain:
-    coordinates = coordinates.tolist()
-    volumes = volumes.tolist()
-    # FIXME:should this be in mesh? (peter row)
-    if fid.smoothing == 'Yes':
-        unique = False
-    else:
-        unique = True
-    if unique:
-        coordinates, volumes, boundary = weed(coordinates, volumes,boundary)
-
-      
-    
-    try:
-        domain = Domain(coordinates, volumes, boundary)
-    except AssertionError, e:
-        fid.close()
-        msg = 'Domain could not be created: %s. ' \
-              'Perhaps use "fail_if_NaN=False and NaN_filler = ..."' % e
-        raise DataDomainError, msg
-
-    if not boundary is None:
-        domain.boundary = boundary
-
-    domain.geo_reference = geo_reference
-
-    domain.starttime = float(starttime) + float(t)
-    domain.time = 0.0
-
-    for quantity in other_quantities:
-        try:
-            NaN = fid.variables[quantity].missing_value
-        except:
-            pass                       # quantity has no missing_value number
-        X = fid.variables[quantity][:]
-        if very_verbose:
-            log.critical('       %s' % str(quantity))
-            log.critical('        NaN = %s' % str(NaN))
-            log.critical('        max(X)')
-            log.critical('       %s' % str(max(X)))
-            log.critical('        max(X)==NaN')
-            log.critical('       %s' % str(max(X)==NaN))
-            log.critical('')
-        if max(X) == NaN or min(X) == NaN:
-            if fail_if_NaN:
-                msg = 'quantity "%s" contains no_data entry' % quantity
-                raise DataMissingValuesError, msg
-            else:
-                data = (X != NaN)
-                X = (X*data) + (data==0)*NaN_filler
-        if unique:
-            X = num.resize(X, (len(X)/3, 3))
-        domain.set_quantity(quantity, X)
-    #
-    for quantity in conserved_quantities:
-        try:
-            NaN = fid.variables[quantity].missing_value
-        except:
-            pass                       # quantity has no missing_value number
-        X = interpolated_quantities[quantity]
-        if very_verbose:
-            log.critical('       %s' % str(quantity))
-            log.critical('        NaN = %s' % str(NaN))
-            log.critical('        max(X)')
-            log.critical('       %s' % str(max(X)))
-            log.critical('        max(X)==NaN')
-            log.critical('       %s' % str(max(X)==NaN))
-            log.critical('')
-        if max(X) == NaN or min(X) == NaN:
-            if fail_if_NaN:
-                msg = 'quantity "%s" contains no_data entry' % quantity
-                raise DataMissingValuesError, msg
-            else:
-                data = (X != NaN)
-                X = (X*data) + (data==0)*NaN_filler
-        if unique:
-            X = num.resize(X, (X.shape[0]/3, 3))
-        domain.set_quantity(quantity, X)
-
-    fid.close()
-
-    return domain
-
-
-##
-# @brief Interpolate a quantity wrt time.
-# @param saved_quantity The quantity to interpolate.
-# @param time_interp (index, ratio)
-# @return A vector of interpolated values.
-def interpolated_quantity(saved_quantity, time_interp):
-    '''Given an index and ratio, interpolate quantity with respect to time.'''
-
-    index, ratio = time_interp
-
-    Q = saved_quantity
-
-    if ratio > 0:
-        q = (1-ratio)*Q[index] + ratio*Q[index+1]
-    else:
-        q = Q[index]
-
-    #Return vector of interpolated values
-    return q
-
-
-##
-# @brief 
-# @parm time 
-# @param t 
-# @return An (index, ration) tuple.
-def get_time_interp(time, t=None):
-    #Finds the ratio and index for time interpolation.
-    #It is borrowed from previous abstract_2d_finite_volumes code.
-    if t is None:
-        t=time[-1]
-        index = -1
-        ratio = 0.
-    else:
-        T = time
-        tau = t
-        index=0
-        msg = 'Time interval derived from file %s [%s:%s]' \
-              % ('FIXMEfilename', T[0], T[-1])
-        msg += ' does not match model time: %s' % tau
-        if tau < time[0]: raise DataTimeError, msg
-        if tau > time[-1]: raise DataTimeError, msg
-        while tau > time[index]: index += 1
-        while tau < time[index]: index -= 1
-        if tau == time[index]:
-            #Protect against case where tau == time[-1] (last time)
-            # - also works in general when tau == time[i]
-            ratio = 0
-        else:
-            #t is now between index and index+1
-            ratio = (tau - time[index])/(time[index+1] - time[index])
-
-    return (index, ratio)
-
-
-##
-# @brief 
-# @param coordinates 
-# @param volumes 
-# @param boundary 
-def weed(coordinates, volumes, boundary=None):
-    if isinstance(coordinates, num.ndarray):
-        coordinates = coordinates.tolist()
-    if isinstance(volumes, num.ndarray):
-        volumes = volumes.tolist()
-
-    unique = False
-    point_dict = {}
-    same_point = {}
-    for i in range(len(coordinates)):
-        point = tuple(coordinates[i])
-        if point_dict.has_key(point):
-            unique = True
-            same_point[i] = point
-            #to change all point i references to point j
-        else:
-            point_dict[point] = i
-            same_point[i] = point
-
-    coordinates = []
-    i = 0
-    for point in point_dict.keys():
-        point = tuple(point)
-        coordinates.append(list(point))
-        point_dict[point] = i
-        i += 1
-
-    for volume in volumes:
-        for i in range(len(volume)):
-            index = volume[i]
-            if index > -1:
-                volume[i] = point_dict[same_point[index]]
-
-    new_boundary = {}
-    if not boundary is None:
-        for segment in boundary.keys():
-            point0 = point_dict[same_point[segment[0]]]
-            point1 = point_dict[same_point[segment[1]]]
-            label = boundary[segment]
-            #FIXME should the bounday attributes be concaterated
-            #('exterior, pond') or replaced ('pond')(peter row)
-
-            if new_boundary.has_key((point0, point1)):
-                new_boundary[(point0,point1)] = new_boundary[(point0, point1)]
-
-            elif new_boundary.has_key((point1, point0)):
-                new_boundary[(point1,point0)] = new_boundary[(point1, point0)]
-            else: new_boundary[(point0, point1)] = label
-
-        boundary = new_boundary
-
-    return coordinates, volumes, boundary
-
-
 ##
 # @brief Read DEM file, decimate data, write new DEM file.
@@ -1164,284 +842,4 @@
 # END URS 2 SWW
 ################################################################################
-
-################################################################################
-# URS UNGRIDDED 2 SWW
-################################################################################
-
-### PRODUCING THE POINTS NEEDED FILE ###
-
-# Ones used for FESA 2007 results
-#LL_LAT = -50.0
-#LL_LONG = 80.0
-#GRID_SPACING = 1.0/60.0
-#LAT_AMOUNT = 4800
-#LONG_AMOUNT = 3600
-
-
-##
-# @brief 
-# @param file_name 
-# @param boundary_polygon 
-# @param zone 
-# @param ll_lat 
-# @param ll_long 
-# @param grid_spacing 
-# @param lat_amount 
-# @param long_amount 
-# @param isSouthernHemisphere 
-# @param export_csv 
-# @param use_cache 
-# @param verbose True if this function is to be verbose.
-# @return 
-def URS_points_needed_to_file(file_name, boundary_polygon, zone,
-                              ll_lat, ll_long,
-                              grid_spacing,
-                              lat_amount, long_amount,
-                              isSouthernHemisphere=True,
-                              export_csv=False, use_cache=False,
-                              verbose=False):
-    """
-    Given the info to replicate the URS grid and a polygon output
-    a file that specifies the cloud of boundary points for URS.
-
-    This creates a .urs file.  This is in the format used by URS;
-    1st line is the number of points,
-    each line after represents a point,in lats and longs.
-
-    Note: The polygon cannot cross zones or hemispheres.
-
-    A work-a-round for different zones or hemispheres is to run this twice,
-    once for each zone, and then combine the output.
-
-    file_name - name of the urs file produced for David.
-    boundary_polygon - a list of points that describes a polygon.
-                      The last point is assumed ot join the first point.
-                      This is in UTM (lat long would be better though)
-
-     This is info about the URS model that needs to be inputted.
-
-    ll_lat - lower left latitude, in decimal degrees
-    ll-long - lower left longitude, in decimal degrees
-    grid_spacing - in deciamal degrees
-    lat_amount - number of latitudes
-    long_amount- number of longs
-
-    Don't add the file extension.  It will be added.
-    """
-
-    geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long,
-                            grid_spacing,
-                            lat_amount, long_amount, isSouthernHemisphere,
-                            use_cache, verbose)
-
-    if not file_name[-4:] == ".urs":
-        file_name += ".urs"
-
-    geo.export_points_file(file_name, isSouthHemisphere=isSouthernHemisphere)
-
-    if export_csv:
-        if file_name[-4:] == ".urs":
-            file_name = file_name[:-4] + ".csv"
-        geo.export_points_file(file_name)
-
-    return geo
-
-
-##
-# @brief 
-# @param boundary_polygon
-# @param zone
-# @param ll_lat
-# @param ll_long
-# @param grid_spacing
-# @param lat_amount
-# @param long_amount
-# @param isSouthHemisphere
-# @param use_cache
-# @param verbose
-def URS_points_needed(boundary_polygon, zone, ll_lat,
-                      ll_long, grid_spacing,
-                      lat_amount, long_amount, isSouthHemisphere=True,
-                      use_cache=False, verbose=False):
-    args = (boundary_polygon,
-            zone, ll_lat,
-            ll_long, grid_spacing,
-            lat_amount, long_amount, isSouthHemisphere)
-    kwargs = {}
-
-    if use_cache is True:
-        try:
-            from anuga.caching import cache
-        except:
-            msg = 'Caching was requested, but caching module' \
-                  'could not be imported'
-            raise msg
-
-        geo = cache(_URS_points_needed,
-                    args, kwargs,
-                    verbose=verbose,
-                    compression=False)
-    else:
-        geo = apply(_URS_points_needed, args, kwargs)
-
-    return geo
-
-
-##
-# @brief 
-# @param boundary_polygon An iterable of points that describe a polygon.
-# @param zone
-# @param ll_lat Lower left latitude, in decimal degrees
-# @param ll_long Lower left longitude, in decimal degrees
-# @param grid_spacing Grid spacing in decimal degrees.
-# @param lat_amount
-# @param long_amount
-# @param isSouthHemisphere
-def _URS_points_needed(boundary_polygon,
-                       zone, ll_lat,
-                       ll_long, grid_spacing,
-                       lat_amount, long_amount,
-                       isSouthHemisphere):
-    """
-    boundary_polygon - a list of points that describes a polygon.
-                      The last point is assumed ot join the first point.
-                      This is in UTM (lat long would b better though)
-
-    ll_lat - lower left latitude, in decimal degrees
-    ll-long - lower left longitude, in decimal degrees
-    grid_spacing - in decimal degrees
-
-    """
-
-    msg = "grid_spacing can not be zero"
-    assert not grid_spacing == 0, msg
-
-    a = boundary_polygon
-
-    # List of segments.  Each segment is two points.
-    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
-
-    # convert the segs to Lat's and longs.
-    # Don't assume the zone of the segments is the same as the lower left
-    # corner of the lat long data!!  They can easily be in different zones
-    lat_long_set = frozenset()
-    for seg in segs:
-        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
-                                        lat_amount, long_amount, zone,
-                                        isSouthHemisphere)
-        lat_long_set |= frozenset(points_lat_long)
-
-    if lat_long_set == frozenset([]):
-        msg = "URS region specified and polygon does not overlap."
-        raise ValueError, msg
-
-    # Warning there is no info in geospatial saying the hemisphere of
-    # these points.  There should be.
-    geo = Geospatial_data(data_points=list(lat_long_set),
-                          points_are_lats_longs=True)
-
-    return geo
-
-
-##
-# @brief Get the points that are needed to interpolate any point a a segment.
-# @param seg Two points in the UTM.
-# @param ll_lat Lower left latitude, in decimal degrees
-# @param ll_long Lower left longitude, in decimal degrees
-# @param grid_spacing 
-# @param lat_amount 
-# @param long_amount 
-# @param zone 
-# @param isSouthHemisphere 
-# @return A list of points.
-def points_needed(seg, ll_lat, ll_long, grid_spacing,
-                  lat_amount, long_amount, zone,
-                  isSouthHemisphere):
-    """
-    seg is two points, in UTM
-    return a list of the points, in lats and longs that are needed to
-    interpolate any point on the segment.
-    """
-
-    from math import sqrt
-
-    geo_reference = Geo_reference(zone=zone)
-    geo = Geospatial_data(seg, geo_reference=geo_reference)
-    seg_lat_long = geo.get_data_points(as_lat_long=True,
-                                       isSouthHemisphere=isSouthHemisphere)
-
-    # 1.415 = 2^0.5, rounded up....
-    sqrt_2_rounded_up = 1.415
-    buffer = sqrt_2_rounded_up * grid_spacing
-
-    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
-    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
-    min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
-    min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
-
-    first_row = (min_long - ll_long) / grid_spacing
-
-    # To round up
-    first_row_long = int(round(first_row + 0.5))
-
-    last_row = (max_long - ll_long) / grid_spacing # round down
-    last_row_long = int(round(last_row))
-
-    first_row = (min_lat - ll_lat) / grid_spacing
-    # To round up
-    first_row_lat = int(round(first_row + 0.5))
-
-    last_row = (max_lat - ll_lat) / grid_spacing # round down
-    last_row_lat = int(round(last_row))
-
-    max_distance = 157147.4112 * grid_spacing
-    points_lat_long = []
-
-    # Create a list of the lat long points to include.
-    for index_lat in range(first_row_lat, last_row_lat + 1):
-        for index_long in range(first_row_long, last_row_long + 1):
-            lat = ll_lat + index_lat*grid_spacing
-            long = ll_long + index_long*grid_spacing
-
-            #filter here to keep good points
-            if keep_point(lat, long, seg, max_distance):
-                points_lat_long.append((lat, long)) #must be hashable
-
-    # Now that we have these points, lets throw ones out that are too far away
-    return points_lat_long
-
-
-##
-# @brief 
-# @param lat
-# @param long
-# @param seg Two points in UTM.
-# @param max_distance
-def keep_point(lat, long, seg, max_distance):
-    """
-    seg is two points, UTM
-    """
-
-    from math import sqrt
-
-    _ , x0, y0 = redfearn(lat, long)
-    x1 = seg[0][0]
-    y1 = seg[0][1]
-    x2 = seg[1][0]
-    y2 = seg[1][1]
-    x2_1 = x2-x1
-    y2_1 = y2-y1
-    try:
-        d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt( \
-            (x2_1)*(x2_1)+(y2_1)*(y2_1))
-    except ZeroDivisionError:
-        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 \
-           and abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
-            return True
-        else:
-            return False
-
-    return d <= max_distance
-
 
 
@@ -1555,89 +953,4 @@
 ################################################################################
 
-##
-# @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
 
 
Index: /trunk/anuga_core/source/anuga/shallow_water/eqf_v2.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/eqf_v2.py	(revision 7768)
+++ /trunk/anuga_core/source/anuga/shallow_water/eqf_v2.py	(revision 7769)
@@ -224,13 +224,13 @@
                     U[i]=U[i]+SIGN*DU[i]
             
-        U1 =U[0]                                                          
-        U2 =U[1]                                                          
-        U3 =U[2]                                                          
-        U11=U[3]                                                          
-        U12=U[4]                                                          
-        U21=U[5]                                                          
-        U22=U[6]                                                          
-        U31=U[7]                                                          
-        U32=U[8]
+        U1  = U[0]                                                          
+        U2  = U[1]                                                          
+        U3  = U[2]                                                          
+        U11 = U[3]                                                          
+        U12 = U[4]                                                          
+        U21 = U[5]                                                          
+        U22 = U[6]                                                          
+        U31 = U[7]                                                          
+        U32 = U[8]
 
         self.U3 = U3
@@ -298,6 +298,6 @@
 
         if CD == 0:
-            #C==============================                                         
-            #C=====   INCLINED FAULT   =====                                         
+            #C==============================     
+            #C=====   INCLINED FAULT   =====     
             #C==============================
             RD2=RD*RD
@@ -311,6 +311,6 @@
             C3= ALP*SD/RD*(XI2*RRD - F1)  
         else:
-            #C==============================                                         
-            #C=====   VERTICAL FAULT   =====                                         
+            #C==============================   
+            #C=====   VERTICAL FAULT   =====     
             #C==============================
             TD=SD/CD                                                          
@@ -345,6 +345,6 @@
         
         if DISL1 != F0:
-            #C======================================                                 
-            #C=====  STRIKE-SLIP CONTRIBUTION  =====                                 
+            #C======================================   
+            #C=====  STRIKE-SLIP CONTRIBUTION  =====  
             #C======================================
             UN=DISL1/PI2                                                      
@@ -362,6 +362,6 @@
 
         if DISL2 != F0:
-            #C===================================                                    
-            #C=====  DIP-SLIP CONTRIBUTION  =====                                    
+            #C===================================  
+            #C=====  DIP-SLIP CONTRIBUTION  =====        
             #C=================================== 
             UN=DISL2/PI2                                                      
@@ -379,5 +379,5 @@
         if DISL3 != F0:
             #C========================================
-            #C=====  TENSILE-FAULT CONTRIBUTION  =====                               
+            #C=====  TENSILE-FAULT CONTRIBUTION  =====   
             #C======================================== 
             UN=DISL3/PI2                                                      
@@ -443,50 +443,50 @@
         PI2 = 6.283185307179586
 
-        D =DEP                                                            
-        P =Y*CD + D*SD                                                    
-        Q =Y*SD - D*CD                                                    
-        S =P*SD + Q*CD                                                    
-        X2=X*X                                                            
-        Y2=Y*Y                                                            
-        XY=X*Y                                                            
-        D2=D*D                                                            
-        R2=X2 + Y2 + D2                                                   
-        R =sqrt(R2)                                                       
-        R3=R *R2
-        R5=R3*R2
-        QR=F3*Q/R5
-        XR =F5*X2/R2
-        YR =F5*Y2/R2
-        XYR=F5*XY/R2
-        DR =F5*D /R2
-        RD =R + D
-        R12=F1/(R*RD*RD)
-        R32=R12*(F2*R + D)/ R2
-        R33=R12*(F3*R + D)/(R2*RD)
-        R53=R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD)
-        R54=R12*(F5*R2 + F4*R*D +    D2)/R3*R12                           
-
-        A1= ALP*Y*(R12-X2*R33)                                            
-        A2= ALP*X*(R12-Y2*R33)                                            
-        A3= ALP*X/R3 - A2                                                 
-        A4=-ALP*XY*R32                                                    
-        A5= ALP*( F1/(R*RD) - X2*R32 )                                    
-        B1= ALP*(-F3*XY*R33      + F3*X2*XY*R54)                          
-        B2= ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54)                          
-        B3= ALP*( F1/R3 - F3*X2/R5) - B2                                  
-        B4=-ALP*F3*XY/R5 - B1                                             
-        C1=-ALP*Y*(R32 - X2*R53)                                          
-        C2=-ALP*X*(R32 - Y2*R53)                                          
-        C3=-ALP*F3*X*D/R5 - C2                                            
-
-        U1 =F0                                                            
-        U2 =F0                                                            
-        U3 =F0                                                            
-        U11=F0                                                            
-        U12=F0                                                            
-        U21=F0                                                            
-        U22=F0                                                            
-        U31=F0                                                            
-        U32=F0
+        D = DEP                                                            
+        P = Y*CD + D*SD                                                    
+        Q = Y*SD - D*CD                                                    
+        S = P*SD + Q*CD                                                    
+        X2 = X*X                                                            
+        Y2 = Y*Y                                                            
+        XY = X*Y                                                            
+        D2 = D*D                                                            
+        R2= X2 + Y2 + D2                                                   
+        R = sqrt(R2)                                                       
+        R3 = R *R2
+        R5 = R3*R2
+        QR = F3*Q/R5
+        XR = F5*X2/R2
+        YR = F5*Y2/R2
+        XYR = F5*XY/R2
+        DR  = F5*D /R2
+        RD = R + D
+        R12 = F1/(R*RD*RD)
+        R32 = R12*(F2*R + D)/ R2
+        R33 = R12*(F3*R + D)/(R2*RD)
+        R53 = R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD)
+        R54 = R12*(F5*R2 + F4*R*D +    D2)/R3*R12                           
+
+        A1 =  ALP*Y*(R12-X2*R33)                                            
+        A2 =  ALP*X*(R12-Y2*R33)                                            
+        A3 =  ALP*X/R3 - A2                                                 
+        A4 = -ALP*XY*R32                                                    
+        A5 =  ALP*( F1/(R*RD) - X2*R32 )                                    
+        B1 =  ALP*(-F3*XY*R33      + F3*X2*XY*R54)                          
+        B2 =  ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54)                          
+        B3 =  ALP*( F1/R3 - F3*X2/R5) - B2                                  
+        B4 = -ALP*F3*XY/R5 - B1                                             
+        C1 = -ALP*Y*(R32 - X2*R53)                                          
+        C2 = -ALP*X*(R32 - Y2*R53)                                          
+        C3 = -ALP*F3*X*D/R5 - C2                                            
+
+        U1 = F0                                                            
+        U2 = F0                                                            
+        U3 = F0                                                            
+        U11= F0                                                            
+        U12= F0                                                            
+        U21= F0                                                            
+        U22= F0                                                            
+        U31= F0                                                            
+        U32= F0
         
         #======================================                                 
@@ -495,54 +495,54 @@
 
         if POT1 != F0:
-            UN=POT1/PI2                                                       
-            QRX=QR*X                                                          
-            FX=F3*X/R5*SD                                                     
-            U1 =U1 - UN*( QRX*X + A1*SD )                                     
-            U2 =U2 - UN*( QRX*Y + A2*SD )                                     
-            U3 =U3 - UN*( QRX*D + A4*SD )
+            UN = POT1/PI2                                                       
+            QRX = QR*X                                                          
+            FX = F3*X/R5*SD                                                     
+            U1 = U1 - UN*( QRX*X + A1*SD )                                     
+            U2 = U2 - UN*( QRX*Y + A2*SD )                                     
+            U3 = U3 - UN*( QRX*D + A4*SD )
               
-            U11=U11- UN*( QRX* (F2-XR)        + B1*SD )                       
-            U12=U12- UN*(-QRX*XYR      + FX*X + B2*SD )                       
-            U21=U21- UN*( QR*Y*(F1-XR)        + B2*SD )                       
-            U22=U22- UN*( QRX *(F1-YR) + FX*Y + B4*SD )                       
-            U31=U31- UN*( QR*D*(F1-XR)        + C1*SD )                       
-            U32=U32- UN*(-QRX*DR*Y     + FX*D + C2*SD )                       
-
-        #======================================                                 
-        #=====    DIP-SLIP CONTRIBUTION   =====                                 
-        #======================================                                 
+            U11 = U11 - UN * ( QRX* (F2-XR)        + B1*SD )  
+            U12 = U12 - UN * (-QRX*XYR      + FX*X + B2*SD )   
+            U21 = U21 - UN * ( QR*Y*(F1-XR)        + B2*SD )   
+            U22 = U22 - UN * ( QRX *(F1-YR) + FX*Y + B4*SD )    
+            U31 = U31 - UN * ( QR*D*(F1-XR)        + C1*SD )  
+            U32 = U32 - UN * (-QRX*DR*Y     + FX*D + C2*SD ) 
+
+        #======================================           
+        #=====    DIP-SLIP CONTRIBUTION   =====          
+        #======================================          
 
         if POT2 != F0:
-            UN=POT2/PI2                                                       
-            SDCD=SD*CD                                                        
-            QRP=QR*P                                                          
-            FS=F3*S/R5                                                        
-            U1 =U1 - UN*( QRP*X - A3*SDCD )                                   
-            U2 =U2 - UN*( QRP*Y - A1*SDCD )                                   
-            U3 =U3 - UN*( QRP*D - A5*SDCD )                                   
-            U11=U11- UN*( QRP*(F1-XR)        - B3*SDCD )                      
-            U12=U12- UN*(-QRP*XYR     + FS*X - B1*SDCD )                      
-            U21=U21- UN*(-QRP*XYR            - B1*SDCD )                      
-            U22=U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD )                      
-            U31=U31- UN*(-QRP*DR*X           - C3*SDCD )                      
-            U32=U32- UN*(-QRP*DR*Y    + FS*D - C1*SDCD )
-            
-        #========================================                               
-        #=====  TENSILE-FAULT CONTRIBUTION  =====                               
-        #========================================                                                                    
+            UN = POT2/PI2                                                       
+            SDCD = SD*CD                                                        
+            QRP = QR*P                                                          
+            FS = F3*S/R5                                                        
+            U1  = U1 - UN*( QRP*X - A3*SDCD )                                   
+            U2  = U2 - UN*( QRP*Y - A1*SDCD )                                   
+            U3  = U3 - UN*( QRP*D - A5*SDCD )                                   
+            U11 = U11- UN*( QRP*(F1-XR)        - B3*SDCD )                      
+            U12 = U12- UN*(-QRP*XYR     + FS*X - B1*SDCD )                      
+            U21 = U21- UN*(-QRP*XYR            - B1*SDCD )                      
+            U22 = U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD )                      
+            U31 = U31- UN*(-QRP*DR*X           - C3*SDCD )                      
+            U32 = U32- UN*(-QRP*DR*Y    + FS*D - C1*SDCD )
+            
+        #========================================    
+        #=====  TENSILE-FAULT CONTRIBUTION  =====    
+        #========================================                             
 
         if POT3 != F0:
-            UN=POT3/PI2                                                       
-            SDSD=SD*SD                                                        
-            QRQ=QR*Q                                                          
-            FQ=F2*QR*SD                                                       
-            U1 =U1 + UN*( QRQ*X - A3*SDSD )                                   
-            U2 =U2 + UN*( QRQ*Y - A1*SDSD )                                   
-            U3 =U3 + UN*( QRQ*D - A5*SDSD )                                   
-            U11=U11+ UN*( QRQ*(F1-XR)        - B3*SDSD )                      
-            U12=U12+ UN*(-QRQ*XYR     + FQ*X - B1*SDSD )                      
-            U21=U21+ UN*(-QRQ*XYR            - B1*SDSD )                      
-            U22=U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD )                      
-            U31=U31+ UN*(-QRQ*DR*X           - C3*SDSD )                      
-            U32=U32+ UN*(-QRQ*DR*Y    + FQ*D - C1*SDSD )
-
+            UN = POT3/PI2                                                       
+            SDSD = SD*SD                                                        
+            QRQ = QR*Q                                                          
+            FQ = F2*QR*SD                                                       
+            U1  = U1 + UN*( QRQ*X - A3*SDSD )                                   
+            U2  = U2 + UN*( QRQ*Y - A1*SDSD )                                   
+            U3  = U3 + UN*( QRQ*D - A5*SDSD )                                   
+            U11 = U11+ UN*( QRQ*(F1-XR)        - B3*SDSD )                      
+            U12 = U12+ UN*(-QRQ*XYR     + FQ*X - B1*SDSD )                      
+            U21 = U21+ UN*(-QRQ*XYR            - B1*SDSD )                      
+            U22 = U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD )                      
+            U31 = U31+ UN*(-QRQ*DR*X           - C3*SDSD )                      
+            U32 = U32+ UN*(-QRQ*DR*Y    + FQ*D - C1*SDSD )
+
Index: /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7768)
+++ /trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7769)
@@ -29,4 +29,5 @@
 from anuga.config import netcdf_float, epsilon, g
 
+
 # import all the boundaries - some are generic, some are shallow water
 from boundaries import Reflective_boundary, \
@@ -44,6 +45,9 @@
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 
-
-class Test_Data_Manager(unittest.TestCase):
+# use helper methods from other unit test
+from anuga.file.test_mux import Test_Mux
+
+
+class Test_Data_Manager(Test_Mux):
     # Class variable
     verbose = False
@@ -1232,138 +1236,4 @@
 
 
-    def test_sww2domain1(self):
-        ################################################
-        #Create a test domain, and evolve and save it.
-        ################################################
-        from mesh_factory import rectangular
-
-        #Create basic mesh
-
-        yiel=0.01
-        points, vertices, boundary = rectangular(10,10)
-
-        #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.geo_reference = Geo_reference(56,11,11)
-        domain.smooth = False
-        domain.store = True
-        domain.set_name('bedslope')
-        domain.default_order=2
-        #Bed-slope and friction
-        domain.set_quantity('elevation', lambda x,y: -x/3)
-        domain.set_quantity('friction', 0.1)
-        # Boundary conditions
-        from math import sin, pi
-        Br = Reflective_boundary(domain)
-        Bt = Transmissive_boundary(domain)
-        Bd = Dirichlet_boundary([0.2,0.,0.])
-        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
-
-        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
-        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
-
-        domain.quantities_to_be_stored['xmomentum'] = 2
-        domain.quantities_to_be_stored['ymomentum'] = 2
-        #Initial condition
-        h = 0.05
-        elevation = domain.quantities['elevation'].vertex_values
-        domain.set_quantity('stage', elevation + h)
-
-        domain.check_integrity()
-        #Evolution
-        #domain.tight_slope_limiters = 1
-        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
-            #domain.write_time()
-            pass
-
-
-        ##########################################
-        #Import the example's file as a new domain
-        ##########################################
-        from data_manager import sww2domain
-        import os
-
-        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
-        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
-        #points, vertices, boundary = rectangular(15,15)
-        #domain2.boundary = boundary
-        ###################
-        ##NOW TEST IT!!!
-        ###################
-
-        os.remove(filename)
-
-        bits = ['vertex_coordinates']
-        for quantity in domain.quantities_to_be_stored:
-            bits.append('get_quantity("%s").get_integral()' % quantity)
-            bits.append('get_quantity("%s").get_values()' % quantity)
-
-        for bit in bits:
-            #print 'testing that domain.'+bit+' has been restored'
-            #print bit
-            #print 'done'
-            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
-
-        ######################################
-        #Now evolve them both, just to be sure
-        ######################################x
-        domain.time = 0.
-        from time import sleep
-
-        final = .1
-        domain.set_quantity('friction', 0.1)
-        domain.store = False
-        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
-
-
-        for t in domain.evolve(yieldstep = yiel, finaltime = final):
-            #domain.write_time()
-            pass
-
-        final = final - (domain2.starttime-domain.starttime)
-        #BUT since domain1 gets time hacked back to 0:
-        final = final + (domain2.starttime-domain.starttime)
-
-        domain2.smooth = False
-        domain2.store = False
-        domain2.default_order=2
-        domain2.set_quantity('friction', 0.1)
-        #Bed-slope and friction
-        # Boundary conditions
-        Bd2=Dirichlet_boundary([0.2,0.,0.])
-        domain2.boundary = domain.boundary
-        #print 'domain2.boundary'
-        #print domain2.boundary
-        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
-        #domain2.set_boundary({'exterior': Bd})
-
-        domain2.check_integrity()
-
-        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
-            #domain2.write_time()
-            pass
-
-        ###################
-        ##NOW TEST IT!!!
-        ##################
-
-        bits = ['vertex_coordinates']
-
-        for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
-            bits.append('get_quantity("%s").get_integral()' %quantity)
-            bits.append('get_quantity("%s").get_values()' %quantity)
-
-        #print bits
-        for bit in bits:
-            #print bit
-            #print eval('domain.'+bit)
-            #print eval('domain2.'+bit)
-            
-            #print eval('domain.'+bit+'-domain2.'+bit)
-            msg = 'Values in the two domains are different for ' + bit
-            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),
-                                rtol=1.e-5, atol=3.e-8), msg
-
-
     def DISABLEDtest_sww2domain2(self):
         ##################################################################
@@ -1413,5 +1283,5 @@
         #Import the file as a new domain
         ##################################
-        from data_manager import sww2domain
+        from data_manager import load_sww_as_domain
         import os
 
@@ -1424,5 +1294,5 @@
         #                     verbose=self.verbose)        
         try:
-            domain2 = sww2domain(filename,
+            domain2 = load_sww_as_domain(filename,
                                  boundary,
                                  fail_if_NaN=True,
@@ -1431,5 +1301,5 @@
             # Now import it, filling NaNs to be -9999
             filler = -9999
-            domain2 = sww2domain(filename,
+            domain2 = load_sww_as_domain(filename,
                                  None,
                                  fail_if_NaN=False,
@@ -1442,5 +1312,5 @@
         # Now import it, filling NaNs to be 0
         filler = -9999
-        domain2 = sww2domain(filename,
+        domain2 = load_sww_as_domain(filename,
                              None,
                              fail_if_NaN=False,
@@ -1551,12 +1421,6 @@
 
 
-        ##########################################
-        #Import the example's file as a new domain
-        ##########################################
-        from data_manager import sww2domain
-        import os
-
         filename = domain.datadir + os.sep + domain.get_name() + '.sww'
-        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
+        domain2 = load_sww_as_domain(filename,None,fail_if_NaN=False,verbose=self.verbose)
         #points, vertices, boundary = rectangular(15,15)
         #domain2.boundary = boundary
@@ -1812,1960 +1676,4 @@
         os.remove(root + '.dem')
         os.remove(root + '_100.dem')      
-        
-    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)
-
-        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)
-        
 
         
@@ -4181,868 +2089,4 @@
     #### END TESTS FOR URS 2 SWW  ###
 
-    #### TESTS URS UNGRIDDED 2 SWW ###
-    def test_URS_points_needed(self):
-        
-        ll_lat = -21.5
-        ll_long = 114.5
-        grid_spacing = 1./60.
-        lat_amount = 30
-        long_amount = 30
-        zone = 50
-
-        boundary_polygon = [[250000,7660000],[280000,7660000],
-                             [280000,7630000],[250000,7630000]]
-        geo=URS_points_needed(boundary_polygon, zone, 
-                              ll_lat, ll_long, grid_spacing, 
-                              lat_amount, long_amount,
-                              verbose=self.verbose)
-        # to test this geo, can info from it be transfered onto the boundary
-        # poly?
-        #Maybe see if I can fit the data to the polygon - have to represent
-        # the poly as points though.
-        #geo.export_points_file("results.txt", as_lat_long=True)
-        results = frozenset(geo.get_data_points(as_lat_long=True))
-        #print 'results',results
-
-        # These are a set of points that have to be in results
-        points = []
-        for i in range(18):
-            lat = -21.0 - 8./60 - grid_spacing * i
-            points.append((lat,degminsec2decimal_degrees(114,35,0))) 
-            points.append((lat,degminsec2decimal_degrees(114,36,0))) 
-            points.append((lat,degminsec2decimal_degrees(114,52,0))) 
-            points.append((lat,degminsec2decimal_degrees(114,53,0)))
-        geo_answer = Geospatial_data(data_points=points,
-                                     points_are_lats_longs=True) 
-        #geo_answer.export_points_file("answer.txt", as_lat_long=True)  
-        answer = frozenset(points)
-        
-        outs = answer.difference(results)
-        #print "outs", outs
-        # This doesn't work.  Though visualising the results shows that
-        # it is correct.
-        #assert answer.issubset(results)
-        # this is why;
-        #point (-21.133333333333333, 114.58333333333333)
-        #result (-21.133333332232368, 114.58333333300342)
-        
-        for point in points:
-            found = False
-            for result in results:
-                if num.allclose(point, result):
-                    found = True
-                    break
-            if not found:
-                assert False
-        
-    
-    def dave_test_URS_points_needed(self):
-        ll_lat = -21.51667
-        ll_long = 114.51667
-        grid_spacing = 2./60.
-        lat_amount = 15
-        long_amount = 15
-
-       
-        boundary_polygon = [[250000,7660000],[280000,7660000],
-                             [280000,7630000],[250000,7630000]]
-        URS_points_needed_to_file('a_test_example',boundary_polygon,
-                                  ll_lat, ll_long, grid_spacing, 
-                                  lat_amount, long_amount,
-                                  verbose=self.verbose)
-        
-    def X_test_URS_points_neededII(self):
-        ll_lat = -21.5
-        ll_long = 114.5
-        grid_spacing = 1./60.
-        lat_amount = 30
-        long_amount = 30
-
-        # change this so lats and longs are inputed, then converted
-        
-        #boundary_polygon = [[7660000,250000],[7660000,280000],
-        #                     [7630000,280000],[7630000,250000]]
-        URS_points_needed(boundary_polygon, ll_lat, ll_long, grid_spacing, 
-                          lat_amount, long_amount,
-                          verbose=self.verbose)
-        
-    def test_URS_points_northern_hemisphere(self):
-               
-        LL_LAT = 8.0
-        LL_LONG = 97.0
-        GRID_SPACING = 2.0/60.0
-        LAT_AMOUNT = 2
-        LONG_AMOUNT = 2
-        ZONE = 47
-
-        # 
-        points = []
-        for i in range(2):
-            for j in range(2):
-                points.append((degminsec2decimal_degrees(8,1+i*2,0),
-                               degminsec2decimal_degrees(97,1+i*2,0)))
-        #print "points", points
-        geo_poly = Geospatial_data(data_points=points,
-                                     points_are_lats_longs=True)
-        poly_lat_long = geo_poly.get_data_points(as_lat_long=False,
-                                       isSouthHemisphere=False)
-        #print "seg_lat_long",  poly_lat_long
-        
-      #   geo=URS_points_needed_to_file('test_example_poly3', poly_lat_long,
-#                                   ZONE,
-#                                   LL_LAT, LL_LONG,
-#                                   GRID_SPACING,
-#                                   LAT_AMOUNT, LONG_AMOUNT,
-#                                   isSouthernHemisphere=False,
-#                                   export_csv=True,
-#                                   verbose=self.verbose)
-        
-        geo=URS_points_needed(poly_lat_long,
-                                  ZONE,
-                                  LL_LAT, LL_LONG,
-                                  GRID_SPACING,
-                                  LAT_AMOUNT, LONG_AMOUNT,
-                                  isSouthHemisphere=False,
-                                  verbose=self.verbose)
-        
-        results = frozenset(geo.get_data_points(as_lat_long=True,
-                                                isSouthHemisphere=False))
-        #print 'results',results
-
-        # These are a set of points that have to be in results
-        points = [] 
-        for i in range(2):
-            for j in range(2):
-                points.append((degminsec2decimal_degrees(8,i*2,0),
-                               degminsec2decimal_degrees(97,i*2,0)))
-        #print "answer points", points
-        answer = frozenset(points)
-        
-        for point in points:
-            found = False
-            for result in results:
-                if num.allclose(point, result):
-                    found = True
-                    break
-            if not found:
-                assert False
-        
-
-    def covered_in_other_tests_test_URS_points_needed_poly1(self):
-        # Values used for FESA 2007 results
-        # domain in southern hemisphere zone 51        
-        LL_LAT = -50.0
-        LL_LONG = 80.0
-        GRID_SPACING = 2.0/60.0
-        LAT_AMOUNT = 4800
-        LONG_AMOUNT = 3600
-        ZONE = 51
-        
-        poly1 = [[296361.89, 8091928.62],
-                 [429495.07,8028278.82],
-                 [447230.56,8000674.05],
-                 [429661.2,7982177.6],
-                 [236945.9,7897453.16],
-                 [183493.44,7942782.27],
-                 [226583.04,8008058.96]]
-
-        URS_points_needed_to_file('test_example_poly2', poly1,
-                                  ZONE,
-                                  LL_LAT, LL_LONG,
-                                  GRID_SPACING,
-                                  LAT_AMOUNT, LONG_AMOUNT,
-                                  verbose=self.verbose)
-        
-
-
-    def covered_in_other_tests_test_URS_points_needed_poly2(self):
-        # Values used for 2004 validation work
-        # domain in northern hemisphere zone 47        
-        LL_LAT = 0.0
-        LL_LONG = 90.0
-        GRID_SPACING = 2.0/60.0
-        LAT_AMOUNT = (15-LL_LAT)/GRID_SPACING
-        LONG_AMOUNT = (100-LL_LONG)/GRID_SPACING 
-        ZONE = 47
-        
-        poly2 = [[419336.424,810100.845],
-                 [342405.0281,711455.8026],
-                 [274649.9152,723352.9603],
-                 [272089.092,972393.0131],
-                 [347633.3754,968551.7784],
-                 [427979.2022,885965.2313],
-                 [427659.0993,875721.9386],
-                 [429259.6138,861317.3083],
-                 [436301.8775,840830.723]]
-        
-        URS_points_needed_to_file('test_example_poly2', poly2,
-                                  ZONE,
-                                  LL_LAT, LL_LONG,
-                                  GRID_SPACING,
-                                  LAT_AMOUNT, LONG_AMOUNT,
-                                  isSouthernHemisphere=False,
-                                  verbose=self.verbose) 
-        
-    #### END TESTS URS UNGRIDDED 2 SWW ###
-    def test_Urs_points(self):
-        time_step_count = 3
-        time_step = 2
-        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
-        base_name, files = self.write_mux(lat_long_points,
-                                          time_step_count, time_step)
-        for file in files:
-            urs = Urs_points(file)
-            assert time_step_count == urs.time_step_count
-            assert time_step == urs.time_step
-
-            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
-                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
-                    assert num.allclose(n, dep[2])
-                        
-            count = 0
-            for slice in urs:
-                count += 1
-                #print slice
-                for lat_lon, quantity in map(None, lat_long_points, slice):
-                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
-                    #print "quantity", quantity
-                    #print "e", e
-                    #print "n", n
-                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
-                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
-                        assert num.allclose(e, quantity)
-                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
-                        assert num.allclose(n, quantity)
-            assert count == time_step_count
-                     
-        self.delete_mux(files)
-
-    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)
-        urs_ungridded2sww(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)
-  
-    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)
-        urs_ungridded2sww(base_name, mean_stage=tide, origin = geo_reference,
-                          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)
-  
-    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)
-        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
-                          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)
-
-        
-    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 '' 
-
-        urs_ungridded2sww(base_name, mean_stage=-240992.0,
-                          hole_points_UTM=[ 292110.784, 7676551.710 ],
-                          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)
-        
-        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 '' 
-
-        urs_ungridded2sww(base_name, mean_stage=-240992.0,
-                          hole_points_UTM=[ 292110.784, 7676551.710 ],
-                          verbose=self.verbose)
-        
-        # 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)
-        
-        urs_ungridded2sww(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)
-        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
-                          mint=101, maxt=500,
-                          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 = [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)
-        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
-                          mint=0, maxt=100000,
-                          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
-        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:
-            urs_ungridded2sww(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:
-            urs_ungridded2sww(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_points_needed(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)
-        urs_ungridded2sww(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)
-        urs_ungridded2sww(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_triangulation(self):
Index: unk/anuga_core/source/anuga/shallow_water/test_file_conversion.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/test_file_conversion.py	(revision 7768)
+++ 	(revision )
@@ -1,2248 +1,0 @@
-
-from shallow_water_domain import Domain
-from file_conversion import ferret2sww
-from anuga.utilities.numerical_tools import ensure_numeric, mean
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
-     import Transmissive_boundary
-from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
-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.config import netcdf_float, epsilon, g
-from Scientific.IO.NetCDF import NetCDFFile
-import sys
-import unittest
-import numpy as num
-import copy
-import os
-
-
-class Test_File_Conversion(unittest.TestCase):
-    """ A suite of tests to test file conversion functions.
-        These tests are quite coarse-grained: converting a file
-        and checking that its headers and some of its contents
-        are correct.
-    """
-    verbose = False
-
-    def set_verbose(self):
-        Test_File_Conversion.verbose = True
-        
-    def setUp(self):
-        import time
-        
-        self.verbose = Test_File_Conversion.verbose
-        # Create basic mesh
-        points, vertices, boundary = rectangular(2, 2)
-
-        # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.default_order = 2
-
-        # Set some field values
-        domain.set_quantity('elevation', lambda x,y: -x)
-        domain.set_quantity('friction', 0.03)
-
-
-        ######################
-        # Boundary conditions
-        B = Transmissive_boundary(domain)
-        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
-
-
-        ######################
-        #Initial condition - with jumps
-        bed = domain.quantities['elevation'].vertex_values
-        stage = num.zeros(bed.shape, num.float)
-
-        h = 0.3
-        for i in range(stage.shape[0]):
-            if i % 2 == 0:
-                stage[i,:] = bed[i,:] + h
-            else:
-                stage[i,:] = bed[i,:]
-
-        domain.set_quantity('stage', stage)
-
-
-        domain.distribute_to_vertices_and_edges()               
-        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
-
-
-        self.domain = domain
-
-        C = domain.get_vertex_coordinates()
-        self.X = C[:,0:6:2].copy()
-        self.Y = C[:,1:6:2].copy()
-
-        self.F = bed
-
-        #Write A testfile (not realistic. Values aren't realistic)
-        self.test_MOST_file = 'most_small'
-
-        longitudes = [150.66667, 150.83334, 151., 151.16667]
-        latitudes = [-34.5, -34.33333, -34.16667, -34]
-
-        long_name = 'LON'
-        lat_name = 'LAT'
-
-        nx = 4
-        ny = 4
-        six = 6
-
-
-        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
-            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
-
-            fid.createDimension(long_name,nx)
-            fid.createVariable(long_name,netcdf_float,(long_name,))
-            fid.variables[long_name].point_spacing='uneven'
-            fid.variables[long_name].units='degrees_east'
-            fid.variables[long_name].assignValue(longitudes)
-
-            fid.createDimension(lat_name,ny)
-            fid.createVariable(lat_name,netcdf_float,(lat_name,))
-            fid.variables[lat_name].point_spacing='uneven'
-            fid.variables[lat_name].units='degrees_north'
-            fid.variables[lat_name].assignValue(latitudes)
-
-            fid.createDimension('TIME',six)
-            fid.createVariable('TIME',netcdf_float,('TIME',))
-            fid.variables['TIME'].point_spacing='uneven'
-            fid.variables['TIME'].units='seconds'
-            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
-
-
-            name = ext[1:3].upper()
-            if name == 'E.': name = 'ELEVATION'
-            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
-            fid.variables[name].units='CENTIMETERS'
-            fid.variables[name].missing_value=-1.e+034
-
-            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
-                                              [-0.1214216, 0, 0, 0],
-                                              [0, 0, 0, 0],
-                                              [0, 0, 0, 0]],
-                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
-                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
-                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
-                                              [0, 2.291054e-005, 2.291054e-005, 0]],
-                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
-                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
-                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
-                                              [0, 0.0001374632, 0.0001374632, 0]],
-                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
-                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
-                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
-                                              [0, 0.0002520159, 0.0002520159, 0]],
-                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
-                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
-                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
-                                              [0, 0.0003665686, 0.0003665686, 0]],
-                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
-                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
-                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
-                                              [0, 0.0004811212, 0.0004811212, 0]]])
-
-
-            fid.close()
-
-
-
-
-    def tearDown(self):
-        import os
-        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
-            #print 'Trying to remove', self.test_MOST_file + ext
-            os.remove(self.test_MOST_file + ext)
-
-    def test_ferret2sww1(self):
-        """Test that georeferencing etc works when converting from
-        ferret format (lat/lon) to sww format (UTM)
-        """
-        from Scientific.IO.NetCDF import NetCDFFile
-        import os, sys
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
-        #
-        # First value (index=0) in small_ha.nc is 0.3400644 cm,
-        # Fourth value (index==3) is -6.50198 cm
-
-
-
-        #Read
-        from anuga.coordinate_transforms.redfearn import redfearn
-        #fid = NetCDFFile(self.test_MOST_file)
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-        first_value = fid.variables['HA'][:][0,0,0]
-        fourth_value = fid.variables['HA'][:][0,0,3]
-        fid.close()
-
-
-        #Call conversion (with zero origin)
-        #ferret2sww('small', verbose=False,
-        #           origin = (56, 0, 0))
-        ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 0, 0))
-
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)
-        #print zone, e, n
-
-        #Read output file 'small.sww'
-        #fid = NetCDFFile('small.sww')
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-
-        #Check that first coordinate is correctly represented
-        assert num.allclose(x[0], e)
-        assert num.allclose(y[0], n)
-
-        #Check first value
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-
-        #print ymomentum
-
-        assert num.allclose(stage[0,0], first_value/100)  #Meters
-
-        #Check fourth value
-        assert num.allclose(stage[0,3], fourth_value/100)  #Meters
-
-        fid.close()
-
-        #Cleanup
-        import os
-        os.remove(self.test_MOST_file + '.sww')
-
-
-
-    def test_ferret2sww_zscale(self):
-        """Test that zscale workse
-        """
-        from Scientific.IO.NetCDF import NetCDFFile
-        import os, sys
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
-        #
-        # First value (index=0) in small_ha.nc is 0.3400644 cm,
-        # Fourth value (index==3) is -6.50198 cm
-
-
-        #Read
-        from anuga.coordinate_transforms.redfearn import redfearn
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-        first_value = fid.variables['HA'][:][0,0,0]
-        fourth_value = fid.variables['HA'][:][0,0,3]
-        fid.close()
-
-        #Call conversion (with no scaling)
-        ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 0, 0))
-
-        #Work out the UTM coordinates for first point
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        #Check values
-        stage_1 = fid.variables['stage'][:]
-        xmomentum_1 = fid.variables['xmomentum'][:]
-        ymomentum_1 = fid.variables['ymomentum'][:]
-
-        assert num.allclose(stage_1[0,0], first_value/100)  #Meters
-        assert num.allclose(stage_1[0,3], fourth_value/100)  #Meters
-
-        fid.close()
-
-        #Call conversion (with scaling)
-        ferret2sww(self.test_MOST_file,
-                   zscale = 5,
-                   verbose=self.verbose,
-                   origin = (56, 0, 0))
-
-        #Work out the UTM coordinates for first point
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        #Check values
-        stage_5 = fid.variables['stage'][:]
-        xmomentum_5 = fid.variables['xmomentum'][:]
-        ymomentum_5 = fid.variables['ymomentum'][:]
-        elevation = fid.variables['elevation'][:]
-
-        assert num.allclose(stage_5[0,0], 5*first_value/100)  #Meters
-        assert num.allclose(stage_5[0,3], 5*fourth_value/100)  #Meters
-
-        assert num.allclose(5*stage_1, stage_5)
-
-        # Momentum will also be changed due to new depth
-
-        depth_1 = stage_1-elevation
-        depth_5 = stage_5-elevation
-
-
-        for i in range(stage_1.shape[0]):
-            for j in range(stage_1.shape[1]):            
-                if depth_1[i,j] > epsilon:
-
-                    scale = depth_5[i,j]/depth_1[i,j]
-                    ref_xmomentum = xmomentum_1[i,j] * scale
-                    ref_ymomentum = ymomentum_1[i,j] * scale
-                    
-                    #print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
-                    
-                    assert num.allclose(xmomentum_5[i,j], ref_xmomentum)
-                    assert num.allclose(ymomentum_5[i,j], ref_ymomentum)
-                    
-        
-
-        fid.close()
-
-
-        #Cleanup
-        import os
-        os.remove(self.test_MOST_file + '.sww')
-
-
-
-    def test_ferret2sww_2(self):
-        """Test that georeferencing etc works when converting from
-        ferret format (lat/lon) to sww format (UTM)
-        """
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
-        #
-        # First value (index=0) in small_ha.nc is 0.3400644 cm,
-        # Fourth value (index==3) is -6.50198 cm
-
-
-        from anuga.coordinate_transforms.redfearn import redfearn
-        #fid = NetCDFFile('small_ha.nc')
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-
-        #Pick a coordinate and a value
-
-        time_index = 1
-        lat_index = 0
-        lon_index = 2
-
-        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
-        test_time = fid.variables['TIME'][:][time_index]
-        test_lat = fid.variables['LAT'][:][lat_index]
-        test_lon = fid.variables['LON'][:][lon_index]
-
-        linear_point_index = lat_index*4 + lon_index
-        fid.close()
-
-        #Call conversion (with zero origin)
-        ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 0, 0))
-
-
-        #Work out the UTM coordinates for test point
-        zone, e, n = redfearn(test_lat, test_lon)
-
-        #Read output file 'small.sww'
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-
-        #Check that test coordinate is correctly represented
-        assert num.allclose(x[linear_point_index], e)
-        assert num.allclose(y[linear_point_index], n)
-
-        #Check test value
-        stage = fid.variables['stage'][:]
-
-        assert num.allclose(stage[time_index, linear_point_index], test_value/100)
-
-        fid.close()
-
-        #Cleanup
-        import os
-        os.remove(self.test_MOST_file + '.sww')
-
-
-    def test_ferret2sww_lat_long(self):
-        # Test that min lat long works
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        
-        #Read
-        from anuga.coordinate_transforms.redfearn import redfearn
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-        first_value = fid.variables['HA'][:][0,0,0]
-        fourth_value = fid.variables['HA'][:][0,0,3]
-        fid.close()
-
-
-        #Call conversion (with zero origin)
-        #ferret2sww('small', verbose=self.verbose,
-        #           origin = (56, 0, 0))
-        ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
-
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)
-        #print zone, e, n
-
-        #Read output file 'small.sww'
-        #fid = NetCDFFile('small.sww')
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        #Check that first coordinate is correctly represented
-        assert 16 == len(x)
-
-        fid.close()
-
-        #Cleanup
-        import os
-        os.remove(self.test_MOST_file + '.sww')
-
-
-    def test_ferret2sww_lat_longII(self):
-        # Test that min lat long works
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        
-        #Read
-        from anuga.coordinate_transforms.redfearn import redfearn
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-        first_value = fid.variables['HA'][:][0,0,0]
-        fourth_value = fid.variables['HA'][:][0,0,3]
-        fid.close()
-
-
-        #Call conversion (with zero origin)
-        #ferret2sww('small', verbose=False,
-        #           origin = (56, 0, 0))
-        ferret2sww(self.test_MOST_file, verbose=False,
-                   origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
-
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)
-        #print zone, e, n
-
-        #Read output file 'small.sww'
-        #fid = NetCDFFile('small.sww')
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        #Check that first coordinate is correctly represented
-        assert 12 == len(x)
-
-        fid.close()
-
-        #Cleanup
-        import os
-        os.remove(self.test_MOST_file + '.sww')
-        
-    def test_ferret2sww3(self):
-        """Elevation included
-        """
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        # ELEVATION = [-1 -2 -3 -4
-        #             -5 -6 -7 -8
-        #              ...
-        #              ...      -16]
-        # where the top left corner is -1m,
-        # and the ll corner is -13.0m
-        #
-        # First value (index=0) in small_ha.nc is 0.3400644 cm,
-        # Fourth value (index==3) is -6.50198 cm
-
-        from anuga.coordinate_transforms.redfearn import redfearn
-        import os
-        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
-        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
-        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
-        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
-
-        h1_list = [150.66667,150.83334,151.]
-        h2_list = [-34.5,-34.33333]
-
-        long_name = 'LON'
-        lat_name = 'LAT'
-        time_name = 'TIME'
-
-        nx = 3
-        ny = 2
-
-        for fid in [fid1,fid2,fid3]:
-            fid.createDimension(long_name,nx)
-            fid.createVariable(long_name,netcdf_float,(long_name,))
-            fid.variables[long_name].point_spacing='uneven'
-            fid.variables[long_name].units='degrees_east'
-            fid.variables[long_name].assignValue(h1_list)
-
-            fid.createDimension(lat_name,ny)
-            fid.createVariable(lat_name,netcdf_float,(lat_name,))
-            fid.variables[lat_name].point_spacing='uneven'
-            fid.variables[lat_name].units='degrees_north'
-            fid.variables[lat_name].assignValue(h2_list)
-
-            fid.createDimension(time_name,2)
-            fid.createVariable(time_name,netcdf_float,(time_name,))
-            fid.variables[time_name].point_spacing='uneven'
-            fid.variables[time_name].units='seconds'
-            fid.variables[time_name].assignValue([0.,1.])
-            #if fid == fid3: break
-
-
-        for fid in [fid4]:
-            fid.createDimension(long_name,nx)
-            fid.createVariable(long_name,netcdf_float,(long_name,))
-            fid.variables[long_name].point_spacing='uneven'
-            fid.variables[long_name].units='degrees_east'
-            fid.variables[long_name].assignValue(h1_list)
-
-            fid.createDimension(lat_name,ny)
-            fid.createVariable(lat_name,netcdf_float,(lat_name,))
-            fid.variables[lat_name].point_spacing='uneven'
-            fid.variables[lat_name].units='degrees_north'
-            fid.variables[lat_name].assignValue(h2_list)
-
-        name = {}
-        name[fid1]='HA'
-        name[fid2]='UA'
-        name[fid3]='VA'
-        name[fid4]='ELEVATION'
-
-        units = {}
-        units[fid1]='cm'
-        units[fid2]='cm/s'
-        units[fid3]='cm/s'
-        units[fid4]='m'
-
-        values = {}
-        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
-        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
-        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
-        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
-
-        for fid in [fid1,fid2,fid3]:
-          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
-          fid.variables[name[fid]].point_spacing='uneven'
-          fid.variables[name[fid]].units=units[fid]
-          fid.variables[name[fid]].assignValue(values[fid])
-          fid.variables[name[fid]].missing_value = -99999999.
-          #if fid == fid3: break
-
-        for fid in [fid4]:
-            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
-            fid.variables[name[fid]].point_spacing='uneven'
-            fid.variables[name[fid]].units=units[fid]
-            fid.variables[name[fid]].assignValue(values[fid])
-            fid.variables[name[fid]].missing_value = -99999999.
-
-
-        fid1.sync(); fid1.close()
-        fid2.sync(); fid2.close()
-        fid3.sync(); fid3.close()
-        fid4.sync(); fid4.close()
-
-        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
-        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
-        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
-
-
-        first_amp = fid1.variables['HA'][:][0,0,0]
-        third_amp = fid1.variables['HA'][:][0,0,2]
-        first_elevation = fid2.variables['ELEVATION'][0,0]
-        third_elevation= fid2.variables['ELEVATION'][:][0,2]
-        first_speed = fid3.variables['VA'][0,0,0]
-        third_speed = fid3.variables['VA'][:][0,0,2]
-
-        fid1.close()
-        fid2.close()
-        fid3.close()
-
-        #Call conversion (with zero origin)
-        ferret2sww('test', verbose=self.verbose,
-                   origin = (56, 0, 0), inverted_bathymetry=False)
-
-        os.remove('test_va.nc')
-        os.remove('test_ua.nc')
-        os.remove('test_ha.nc')
-        os.remove('test_e.nc')
-
-        #Read output file 'test.sww'
-        fid = NetCDFFile('test.sww')
-
-
-        #Check first value
-        elevation = fid.variables['elevation'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-
-        #print ymomentum
-        first_height = first_amp/100 - first_elevation
-        third_height = third_amp/100 - third_elevation
-        first_momentum=first_speed*first_height/100
-        third_momentum=third_speed*third_height/100
-
-        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
-        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
-
-        fid.close()
-
-        #Cleanup
-        os.remove('test.sww')
-
-
-
-    def test_ferret2sww4(self):
-        """Like previous but with augmented variable names as
-        in files produced by ferret as opposed to MOST
-        """
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        # ELEVATION = [-1 -2 -3 -4
-        #             -5 -6 -7 -8
-        #              ...
-        #              ...      -16]
-        # where the top left corner is -1m,
-        # and the ll corner is -13.0m
-        #
-        # First value (index=0) in small_ha.nc is 0.3400644 cm,
-        # Fourth value (index==3) is -6.50198 cm
-
-        from anuga.coordinate_transforms.redfearn import redfearn
-        import os
-        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
-        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
-        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
-        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
-
-        h1_list = [150.66667,150.83334,151.]
-        h2_list = [-34.5,-34.33333]
-
-#        long_name = 'LON961_1261'
-#        lat_name = 'LAT481_841'
-#        time_name = 'TIME1'
-
-        long_name = 'LON'
-        lat_name = 'LAT'
-        time_name = 'TIME'
-
-        nx = 3
-        ny = 2
-
-        for fid in [fid1,fid2,fid3]:
-            fid.createDimension(long_name,nx)
-            fid.createVariable(long_name,netcdf_float,(long_name,))
-            fid.variables[long_name].point_spacing='uneven'
-            fid.variables[long_name].units='degrees_east'
-            fid.variables[long_name].assignValue(h1_list)
-
-            fid.createDimension(lat_name,ny)
-            fid.createVariable(lat_name,netcdf_float,(lat_name,))
-            fid.variables[lat_name].point_spacing='uneven'
-            fid.variables[lat_name].units='degrees_north'
-            fid.variables[lat_name].assignValue(h2_list)
-
-            fid.createDimension(time_name,2)
-            fid.createVariable(time_name,netcdf_float,(time_name,))
-            fid.variables[time_name].point_spacing='uneven'
-            fid.variables[time_name].units='seconds'
-            fid.variables[time_name].assignValue([0.,1.])
-            #if fid == fid3: break
-
-
-        for fid in [fid4]:
-            fid.createDimension(long_name,nx)
-            fid.createVariable(long_name,netcdf_float,(long_name,))
-            fid.variables[long_name].point_spacing='uneven'
-            fid.variables[long_name].units='degrees_east'
-            fid.variables[long_name].assignValue(h1_list)
-
-            fid.createDimension(lat_name,ny)
-            fid.createVariable(lat_name,netcdf_float,(lat_name,))
-            fid.variables[lat_name].point_spacing='uneven'
-            fid.variables[lat_name].units='degrees_north'
-            fid.variables[lat_name].assignValue(h2_list)
-
-        name = {}
-        name[fid1]='HA'
-        name[fid2]='UA'
-        name[fid3]='VA'
-        name[fid4]='ELEVATION'
-
-        units = {}
-        units[fid1]='cm'
-        units[fid2]='cm/s'
-        units[fid3]='cm/s'
-        units[fid4]='m'
-
-        values = {}
-        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
-        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
-        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
-        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
-
-        for fid in [fid1,fid2,fid3]:
-          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
-          fid.variables[name[fid]].point_spacing='uneven'
-          fid.variables[name[fid]].units=units[fid]
-          fid.variables[name[fid]].assignValue(values[fid])
-          fid.variables[name[fid]].missing_value = -99999999.
-          #if fid == fid3: break
-
-        for fid in [fid4]:
-            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
-            fid.variables[name[fid]].point_spacing='uneven'
-            fid.variables[name[fid]].units=units[fid]
-            fid.variables[name[fid]].assignValue(values[fid])
-            fid.variables[name[fid]].missing_value = -99999999.
-
-
-        fid1.sync(); fid1.close()
-        fid2.sync(); fid2.close()
-        fid3.sync(); fid3.close()
-        fid4.sync(); fid4.close()
-
-        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
-        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
-        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
-
-
-        first_amp = fid1.variables['HA'][:][0,0,0]
-        third_amp = fid1.variables['HA'][:][0,0,2]
-        first_elevation = fid2.variables['ELEVATION'][0,0]
-        third_elevation= fid2.variables['ELEVATION'][:][0,2]
-        first_speed = fid3.variables['VA'][0,0,0]
-        third_speed = fid3.variables['VA'][:][0,0,2]
-
-        fid1.close()
-        fid2.close()
-        fid3.close()
-
-        #Call conversion (with zero origin)
-        ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
-                   , inverted_bathymetry=False)
-
-        os.remove('test_va.nc')
-        os.remove('test_ua.nc')
-        os.remove('test_ha.nc')
-        os.remove('test_e.nc')
-
-        #Read output file 'test.sww'
-        fid = NetCDFFile('test.sww')
-
-
-        #Check first value
-        elevation = fid.variables['elevation'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-
-        #print ymomentum
-        first_height = first_amp/100 - first_elevation
-        third_height = third_amp/100 - third_elevation
-        first_momentum=first_speed*first_height/100
-        third_momentum=third_speed*third_height/100
-
-        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
-        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
-
-        fid.close()
-
-        #Cleanup
-        os.remove('test.sww')
-
-
-
-
-    def test_ferret2sww_nz_origin(self):
-        from Scientific.IO.NetCDF import NetCDFFile
-        from anuga.coordinate_transforms.redfearn import redfearn
-
-        #Call conversion (with nonzero origin)
-        ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 100000, 200000))
-
-
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)
-
-        #Read output file 'small.sww'
-        #fid = NetCDFFile('small.sww', netcdf_mode_r)
-        fid = NetCDFFile(self.test_MOST_file + '.sww')
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-
-        #Check that first coordinate is correctly represented
-        assert num.allclose(x[0], e-100000)
-        assert num.allclose(y[0], n-200000)
-
-        fid.close()
-
-        #Cleanup
-        os.remove(self.test_MOST_file + '.sww')
-
-
-    def test_ferret2sww_lat_longII(self):
-        # Test that min lat long works
-
-        #The test file has
-        # LON = 150.66667, 150.83334, 151, 151.16667
-        # LAT = -34.5, -34.33333, -34.16667, -34 ;
-        
-        #Read
-        from anuga.coordinate_transforms.redfearn import redfearn
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
-        first_value = fid.variables['HA'][:][0,0,0]
-        fourth_value = fid.variables['HA'][:][0,0,3]
-        fid.close()
-
-
-        #Call conversion (with zero origin)
-        #ferret2sww('small', verbose=self.verbose,
-        #           origin = (56, 0, 0))
-        try:
-            ferret2sww(self.test_MOST_file, verbose=self.verbose,
-                   origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
-        except AssertionError:
-            pass
-        else:
-            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
-
-    def test_sww_extent(self):
-        """Not a test, rather a look at the sww format
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        self.domain.set_name('datatest' + str(id(self)))
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.reduction = mean
-        self.domain.set_datadir('.')
-        #self.domain.tight_slope_limiters = 1        
-
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        self.domain.time = 2.
-
-        #Modify stage at second timestep
-        stage = self.domain.quantities['stage'].vertex_values
-        self.domain.set_quantity('stage', stage/2)
-
-        sww.store_timestep()
-
-        file_and_extension_name = self.domain.get_name() + ".sww"
-        #print "file_and_extension_name",file_and_extension_name
-        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
-               extent_sww(file_and_extension_name )
-
-        assert num.allclose(xmin, 0.0)
-        assert num.allclose(xmax, 1.0)
-        assert num.allclose(ymin, 0.0)
-        assert num.allclose(ymax, 1.0)
-
-        # FIXME (Ole): Revisit these numbers
-        #assert num.allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
-        #assert num.allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
-
-
-        #Cleanup
-        os.remove(sww.filename)
-
-
-
-#-------------------------------------------------------------
-
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_File_Conversion, 'test_sww')
-    
-    # FIXME(Ole): When Ross has implemented logging, we can 
-    # probably get rid of all this:
-    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
-        Test_File_Conversion.verbose=True
-        saveout = sys.stdout   
-        filename = ".temp_verbose"
-        fid = open(filename, 'w')
-        sys.stdout = fid
-    else:
-        pass
-    runner = unittest.TextTestRunner() #verbosity=2)
-    runner.run(suite)
-
-    # Cleaning up
-    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
-        sys.stdout = saveout 
-        fid.close() 
-        os.remove(filename)
-
-
-    def test_asc_csiro2sww(self):
-        import tempfile
-
-        bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
-        #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( )
-        elevation_dir =  tempfile.mkdtemp()
-        #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
-
-        fid = open(bath_dir_filename, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    0.000    3000.0
-     0.000     9000.000     0.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000
-""")
-        fid.close()
-
-        ucur_dir =  tempfile.mkdtemp()
-        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
-        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-
-        fid = open(ucur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(ucur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        vcur_dir =  tempfile.mkdtemp()
-        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
-        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-
-        fid = open(vcur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(vcur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        sww_file = 'a_test.sww'
-        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
-                      verbose=self.verbose)
-
-        # check the sww file
-
-        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref
-        x_ref = geo_ref.get_xllcorner()
-        y_ref = geo_ref.get_yllcorner()
-        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
-        assert num.allclose(x_ref, 587798.418) # (-38, 148)
-        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
-
-        #Zone:   55
-        #Easting:  588095.674  Northing: 5821451.722
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
-        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
-
-        #Zone:   55
-        #Easting:  632145.632  Northing: 5820863.269
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
-        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
-
-        #Zone:   55
-        #Easting:  609748.788  Northing: 5793447.860
-        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
-        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
-
-        assert num.allclose(z[0],9000.0 )
-        assert num.allclose(stage[0][1],0.0 )
-
-        #(4000+1000)*60
-        assert num.allclose(xmomentum[1][1],300000.0 )
-
-
-        fid.close()
-
-        #tidy up
-        os.remove(bath_dir_filename)
-        os.rmdir(bath_dir)
-
-        os.remove(elevation_dir_filename1)
-        os.remove(elevation_dir_filename2)
-        os.rmdir(elevation_dir)
-
-        os.remove(ucur_dir_filename1)
-        os.remove(ucur_dir_filename2)
-        os.rmdir(ucur_dir)
-
-        os.remove(vcur_dir_filename1)
-        os.remove(vcur_dir_filename2)
-        os.rmdir(vcur_dir)
-
-
-        # remove sww file
-        os.remove(sww_file)
-
-    def test_asc_csiro2sww2(self):
-        import tempfile
-
-        bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
-        #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( )
-        elevation_dir =  tempfile.mkdtemp()
-        #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
-
-        fid = open(bath_dir_filename, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    0.000    3000.0
-     0.000     -9999.000     -9999.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000
-""")
-        fid.close()
-
-        ucur_dir =  tempfile.mkdtemp()
-        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
-        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-
-        fid = open(ucur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(ucur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        vcur_dir =  tempfile.mkdtemp()
-        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
-        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-
-        fid = open(vcur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(vcur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        try:
-            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
-                          vcur_dir, sww_file,
-                      verbose=self.verbose)
-        except:
-            #tidy up
-            os.remove(bath_dir_filename)
-            os.rmdir(bath_dir)
-
-            os.remove(elevation_dir_filename1)
-            os.remove(elevation_dir_filename2)
-            os.rmdir(elevation_dir)
-
-            os.remove(ucur_dir_filename1)
-            os.remove(ucur_dir_filename2)
-            os.rmdir(ucur_dir)
-
-            os.remove(vcur_dir_filename1)
-            os.remove(vcur_dir_filename2)
-            os.rmdir(vcur_dir)
-        else:
-            #tidy up
-            os.remove(bath_dir_filename)
-            os.rmdir(bath_dir)
-
-            os.remove(elevation_dir_filename1)
-            os.remove(elevation_dir_filename2)
-            os.rmdir(elevation_dir)
-            raise 'Should raise exception'
-
-            os.remove(ucur_dir_filename1)
-            os.remove(ucur_dir_filename2)
-            os.rmdir(ucur_dir)
-
-            os.remove(vcur_dir_filename1)
-            os.remove(vcur_dir_filename2)
-            os.rmdir(vcur_dir)
-
-
-
-    def test_asc_csiro2sww3(self):
-        import tempfile
-
-        bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
-        #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( )
-        elevation_dir =  tempfile.mkdtemp()
-        #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
-
-        fid = open(bath_dir_filename, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    0.000    3000.0
-     0.000     -9999.000     -9999.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000
-""")
-        fid.close()
-
-        ucur_dir =  tempfile.mkdtemp()
-        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
-        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-
-        fid = open(ucur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(ucur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        vcur_dir =  tempfile.mkdtemp()
-        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
-        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-
-        fid = open(vcur_dir_filename1, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-        fid = open(vcur_dir_filename2, 'w')
-        fid.write(""" ncols             3
- nrows             2
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-    90.000    60.000    30.0
-    10.000    10.000    10.000
-""")
-        fid.close()
-
-        sww_file = 'a_test.sww'
-        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
-                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
-                      mean_stage = 100,
-                      verbose=self.verbose)
-
-        # check the sww file
-
-        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref
-        x_ref = geo_ref.get_xllcorner()
-        y_ref = geo_ref.get_yllcorner()
-        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
-        assert num.allclose(x_ref, 587798.418) # (-38, 148)
-        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
-
-        #Zone:   55
-        #Easting:  588095.674  Northing: 5821451.722
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
-        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
-
-        #Zone:   55
-        #Easting:  632145.632  Northing: 5820863.269
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
-        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
-
-        #Zone:   55
-        #Easting:  609748.788  Northing: 5793447.860
-        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
-        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
-
-        assert num.allclose(z[0],9000.0 )
-        assert num.allclose(stage[0][4],100.0 )
-        assert num.allclose(stage[0][5],100.0 )
-
-        #(100.0 - 9000)*10
-        assert num.allclose(xmomentum[0][4], -89000.0 )
-
-        #(100.0 - -1000.000)*10
-        assert num.allclose(xmomentum[0][5], 11000.0 )
-
-        fid.close()
-
-        #tidy up
-        os.remove(bath_dir_filename)
-        os.rmdir(bath_dir)
-
-        os.remove(elevation_dir_filename1)
-        os.remove(elevation_dir_filename2)
-        os.rmdir(elevation_dir)
-
-        os.remove(ucur_dir_filename1)
-        os.remove(ucur_dir_filename2)
-        os.rmdir(ucur_dir)
-
-        os.remove(vcur_dir_filename1)
-        os.remove(vcur_dir_filename2)
-        os.rmdir(vcur_dir)
-
-        # remove sww file
-        os.remove(sww_file)
-
-
-    def test_asc_csiro2sww4(self):
-        """
-        Test specifying the extent
-        """
-
-        import tempfile
-
-        bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
-        #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( )
-        elevation_dir =  tempfile.mkdtemp()
-        #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
-
-        fid = open(bath_dir_filename, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -9000.000    -1000.000   -3000.0 -2000.000
-   -1000.000    9000.000  -1000.000 -3000.000
-   -4000.000    6000.000   2000.000 -5000.000
-   -9000.000    -1000.000   -3000.0 -2000.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename1, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -900.000    -100.000   -300.0 -200.000
-   -100.000    900.000  -100.000 -300.000
-   -400.000    600.000   200.000 -500.000
-   -900.000    -100.000   -300.0 -200.000
-""")
-        fid.close()
-
-        fid = open(elevation_dir_filename2, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -990.000    -110.000   -330.0 -220.000
-   -110.000    990.000  -110.000 -330.000
-   -440.000    660.000   220.000 -550.000
-   -990.000    -110.000   -330.0 -220.000
-""")
-        fid.close()
-
-        ucur_dir =  tempfile.mkdtemp()
-        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
-        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-
-        fid = open(ucur_dir_filename1, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -90.000    -10.000   -30.0 -20.000
-   -10.000    90.000  -10.000 -30.000
-   -40.000    60.000   20.000 -50.000
-   -90.000    -10.000   -30.0 -20.000
-""")
-        fid.close()
-        fid = open(ucur_dir_filename2, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -90.000    -10.000   -30.0 -20.000
-   -10.000    99.000  -11.000 -30.000
-   -40.000    66.000   22.000 -50.000
-   -90.000    -10.000   -30.0 -20.000
-""")
-        fid.close()
-
-        vcur_dir =  tempfile.mkdtemp()
-        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
-        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-
-        fid = open(vcur_dir_filename1, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -90.000    -10.000   -30.0 -20.000
-   -10.000    80.000  -20.000 -30.000
-   -40.000    50.000   10.000 -50.000
-   -90.000    -10.000   -30.0 -20.000
-""")
-        fid.close()
-        fid = open(vcur_dir_filename2, 'w')
-        fid.write(""" ncols             4
- nrows             4
- xllcorner    148.00000
- yllcorner    -38.00000
- cellsize       0.25
- nodata_value   -9999.0
-   -90.000    -10.000   -30.0 -20.000
-   -10.000    88.000  -22.000 -30.000
-   -40.000    55.000   11.000 -50.000
-   -90.000    -10.000   -30.0 -20.000
-""")
-        fid.close()
-
-        sww_file = tempfile.mktemp(".sww")
-        #sww_file = 'a_test.sww'
-        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
-                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
-                      mean_stage = 100,
-                       minlat = -37.6, maxlat = -37.6,
-                  minlon = 148.3, maxlon = 148.3,
-                      verbose=self.verbose
-                      #,verbose = True
-                      )
-
-        # check the sww file
-
-        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-        geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref
-        x_ref = geo_ref.get_xllcorner()
-        y_ref = geo_ref.get_yllcorner()
-        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
-
-        assert num.allclose(fid.starttime, 0.0) # (-37.45, 148.25)
-        assert num.allclose(x_ref, 610120.388) # (-37.45, 148.25)
-        assert num.allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
-
-        #Easting:  632145.632  Northing: 5820863.269
-        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
-
-        #print "x",x
-        #print "y",y
-        self.failUnless(len(x) == 4,'failed') # 2*2
-        self.failUnless(len(x) == 4,'failed') # 2*2
-
-        #Zone:   55
-        #Easting:  632145.632  Northing: 5820863.269
-        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
-        # magic number - y is close enough for me.
-        assert num.allclose(x[3], 632145.63 - x_ref)
-        assert num.allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
-
-        assert num.allclose(z[0],9000.0 ) #z is elevation info
-        #print "z",z
-        # 2 time steps, 4 points
-        self.failUnless(xmomentum.shape == (2,4), 'failed')
-        self.failUnless(ymomentum.shape == (2,4), 'failed')
-
-        #(100.0 - -1000.000)*10
-        #assert num.allclose(xmomentum[0][5], 11000.0 )
-
-        fid.close()
-
-        # is the sww file readable?
-        #Lets see if we can convert it to a dem!
-        # if you uncomment, remember to delete the file
-        #print "sww_file",sww_file
-        #dem_file = tempfile.mktemp(".dem")
-        domain = sww2domain(sww_file) ###, dem_file)
-        domain.check_integrity()
-
-        #tidy up
-        os.remove(bath_dir_filename)
-        os.rmdir(bath_dir)
-
-        os.remove(elevation_dir_filename1)
-        os.remove(elevation_dir_filename2)
-        os.rmdir(elevation_dir)
-
-        os.remove(ucur_dir_filename1)
-        os.remove(ucur_dir_filename2)
-        os.rmdir(ucur_dir)
-
-        os.remove(vcur_dir_filename1)
-        os.remove(vcur_dir_filename2)
-        os.rmdir(vcur_dir)
-
-        # remove sww file
-        os.remove(sww_file)
-
-
-
-    def test_get_min_max_indexes(self):
-        latitudes = [3,2,1,0]
-        longitudes = [0,10,20,30]
-
-        # k - lat
-        # l - lon
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            -10,4,-10,31)
-
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-        self.failUnless(latitudes == latitudes_new and \
-                        longitudes == longitudes_news,
-                         'failed')
-
-        ## 2nd test
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            0.5,2.5,5,25)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes == latitudes_new and \
-                        longitudes == longitudes_news,
-                         'failed')
-
-        ## 3rd test
-        kmin, kmax, lmin, lmax = get_min_max_indices(\
-            latitudes,
-            longitudes,
-            1.1,1.9,12,17)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [2, 1] and \
-                        longitudes_news == [10, 20],
-                         'failed')
-
-
-        ## 4th test
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-                                                      -0.1,1.9,-2,17)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [2, 1, 0] and \
-                        longitudes_news == [0, 10, 20],
-                         'failed')
-        ## 5th test
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            0.1,1.9,2,17)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [2, 1, 0] and \
-                        longitudes_news == [0, 10, 20],
-                         'failed')
-
-        ## 6th test
-
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            1.5,4,18,32)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [3, 2, 1] and \
-                        longitudes_news == [10, 20, 30],
-                         'failed')
-
-
-        ## 7th test
-        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            1.5,1.5,15,15)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        m2d = m2d[kmin:kmax,lmin:lmax]
-        #print "m2d", m2d
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(num.alltrue(latitudes_new == [2, 1]) and 
-                        num.alltrue(longitudes_news == [10, 20]),
-                        'failed')
-
-        self.failUnless(num.alltrue(m2d == [[5,6],[9,10]]), 'failed')
-
-    def test_get_min_max_indexes_lat_ascending(self):
-        latitudes = [0,1,2,3]
-        longitudes = [0,10,20,30]
-
-        # k - lat
-        # l - lon
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            -10,4,-10,31)
-
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-        self.failUnless(latitudes == latitudes_new and \
-                        longitudes == longitudes_news,
-                         'failed')
-
-        ## 3rd test
-        kmin, kmax, lmin, lmax = get_min_max_indices(\
-            latitudes,
-            longitudes,
-            1.1,1.9,12,17)
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [1, 2] and \
-                        longitudes_news == [10, 20],
-                         'failed')
-
-    def test_get_min_max_indexes2(self):
-        latitudes = [-30,-35,-40,-45]
-        longitudes = [148,149,150,151]
-
-        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
-
-        # k - lat
-        # l - lon
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            -37,-27,147,149.5)
-
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        #print "m2d", m2d
-        #print "latitudes", latitudes
-        #print "longitudes",longitudes
-        #print "latitudes[kmax]", latitudes[kmax]
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_new = longitudes[lmin:lmax]
-        m2d = m2d[kmin:kmax,lmin:lmax]
-        #print "m2d", m2d
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_new",longitudes_new
-
-        self.failUnless(latitudes_new == [-30, -35, -40] and
-                        longitudes_new == [148, 149,150],
-                        'failed')
-        self.failUnless(num.alltrue(m2d == [[0,1,2],[4,5,6],[8,9,10]]),
-                        'failed')
-
-    def test_get_min_max_indexes3(self):
-        latitudes = [-30,-35,-40,-45,-50,-55,-60]
-        longitudes = [148,149,150,151]
-
-        # k - lat
-        # l - lon
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes,
-            -43,-37,148.5,149.5)
-
-
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        #print "latitudes", latitudes
-        #print "longitudes",longitudes
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == [-35, -40, -45] and
-                        longitudes_news == [148, 149,150],
-                         'failed')
-
-    def test_get_min_max_indexes4(self):
-        latitudes = [-30,-35,-40,-45,-50,-55,-60]
-        longitudes = [148,149,150,151]
-
-        # k - lat
-        # l - lon
-        kmin, kmax, lmin, lmax = get_min_max_indices(
-            latitudes,longitudes)
-
-        #print "kmin",kmin;print "kmax",kmax
-        #print "lmin",lmin;print "lmax",lmax
-        #print "latitudes", latitudes
-        #print "longitudes",longitudes
-        latitudes_new = latitudes[kmin:kmax]
-        longitudes_news = longitudes[lmin:lmax]
-        #print "latitudes_new", latitudes_new
-        #print "longitudes_news",longitudes_news
-
-        self.failUnless(latitudes_new == latitudes  and \
-                        longitudes_news == longitudes,
-                         'failed')
-
-    def test_tsh2sww(self):
-        import os
-        import tempfile
-
-        tsh_file = tempfile.mktemp(".tsh")
-        file = open(tsh_file,"w")
-        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
-0 0.0 0.0 0.0 0.0 0.01 \n \
-1 1.0 0.0 10.0 10.0 0.02  \n \
-2 0.0 1.0 0.0 10.0 0.03  \n \
-3 0.5 0.25 8.0 12.0 0.04  \n \
-# Vert att title  \n \
-elevation  \n \
-stage  \n \
-friction  \n \
-2 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
-0 0 3 2 -1  -1  1 dsg\n\
-1 0 1 3 -1  0 -1   ole nielsen\n\
-4 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
-0 1 0 2 \n\
-1 0 2 3 \n\
-2 2 3 \n\
-3 3 1 1 \n\
-3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
-0 216.0 -86.0 \n \
-1 160.0 -167.0 \n \
-2 114.0 -91.0 \n \
-3 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
-0 0 1 0 \n \
-1 1 2 0 \n \
-2 2 0 0 \n \
-0 # <x> <y> ...Mesh Holes... \n \
-0 # <x> <y> <attribute>...Mesh Regions... \n \
-0 # <x> <y> <attribute>...Mesh Regions, area... \n\
-#Geo reference \n \
-56 \n \
-140 \n \
-120 \n")
-        file.close()
-
-        #sww_file = tempfile.mktemp(".sww")
-        #print "sww_file",sww_file
-        #print "sww_file",tsh_file
-        tsh2sww(tsh_file,
-                verbose=self.verbose)
-
-        os.remove(tsh_file)
-        os.remove(tsh_file[:-4] + '.sww')
-        
-          
-    def test_urs2sww_test_fail(self):
-        points_num = -100
-        time_step_count = 45
-        time_step = -7
-        file_handle, base_name = tempfile.mkstemp("")        
-        os.close(file_handle)
-        os.remove(base_name)
-        
-        files = []
-        quantities = ['HA','UA','VA']
-        
-        mux_names = [WAVEHEIGHT_MUX_LABEL,
-                     EAST_VELOCITY_LABEL,
-                     NORTH_VELOCITY_LABEL]
-        for i,q in enumerate(quantities): 
-            #Write C files
-            columns = 3 # long, lat , depth
-            file = base_name + mux_names[i]
-            f = open(file, 'wb')
-            files.append(file)
-            f.write(pack('i',points_num))
-            f.write(pack('i',time_step_count))
-            f.write(pack('f',time_step))
-
-            f.close()   
-        tide = 1
-        try:
-            urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
-                      verbose=self.verbose)        
-        except ANUGAError:
-            pass
-        else:
-            self.delete_mux(files)
-            msg = 'Should have raised exception'
-            raise msg
-        sww_file = base_name + '.sww'
-        self.delete_mux(files)
-        
-    def test_urs2sww_test_fail2(self):
-        base_name = 'Harry-high-pants'
-        try:
-            urs2sww(base_name)        
-        except IOError:
-            pass
-        else:
-            self.delete_mux(files)
-            msg = 'Should have raised exception'
-            raise msg
-           
-    def test_urs2sww(self):
-        tide = 1
-        base_name, files = self.create_mux()
-        urs2sww(base_name
-                #, origin=(0,0,0)
-                , mean_stage=tide
-                , remove_nc_files=True,
-                      verbose=self.verbose
-                )
-        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)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        geo_reference = Geo_reference(NetCDFObject=fid)
-
-        
-        #Check that first coordinate is correctly represented       
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)
-       
-        assert num.allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
-
-        # Make x and y absolute
-        points = geo_reference.get_absolute(map(None, x, y))
-        points = ensure_numeric(points)
-        x = points[:,0]
-        y = points[:,1]
-        
-        #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) since elevation is inverted in mux files
-        #momentum = velocity_va *(stage+elevation)
-        # = e*(e+tide+n) based on how I'm writing these files
-        answer_y = e*(e+tide+n) * -1 # talking into account mux file format
-        actual_y = ymomentum[0,0]
-        #print "answer_y",answer_y
-        #print "actual_y",actual_y 
-        assert num.allclose(answer_y, actual_y)  #Meters
-        
-        assert num.allclose(answer_x, actual_x)  #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_urs2sww_momentum(self):
-        tide = 1
-        time_step_count = 3
-        time_step = 2
-        #lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
-        # This is gridded
-        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
-        depth=20
-        ha=2
-        ua=5
-        va=-10 #-ve added to take into account mux file format where south
-               # is positive.
-        base_name, files = self.write_mux(lat_long_points,
-                                          time_step_count, time_step,
-                                          depth=depth,
-                                          ha=ha,
-                                          ua=ua,
-                                          va=va)
-        # write_mux(self,lat_long_points, time_step_count, time_step,
-        #          depth=None, ha=None, ua=None, va=None
-        urs2sww(base_name
-                #, origin=(0,0,0)
-                , mean_stage=tide
-                , remove_nc_files=True,
-                      verbose=self.verbose
-                )
-        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)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        geo_reference = Geo_reference(NetCDFObject=fid)
-        
-        #Check first value
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-        elevation = fid.variables['elevation'][:]
-        #assert allclose(stage[0,0], e + tide)  #Meters
-        #print "xmomentum", xmomentum
-        #print "ymomentum", ymomentum
-        #Check the momentums - ua
-        #momentum = velocity*water height
-        #water height = mux_depth + mux_height +tide
-        #water height = mux_depth + mux_height +tide
-        #momentum = velocity*(mux_depth + mux_height +tide)
-        #
-        
-        answer = 115
-        actual = xmomentum[0,0]
-        assert num.allclose(answer, actual)  #Meters^2/ sec
-        answer = 230
-        actual = ymomentum[0,0]
-        #print "answer",answer
-        #print "actual",actual 
-        assert num.allclose(answer, actual)  #Meters^2/ sec
-
-        # check the stage values, first time step.
-        # These arrays are equal since the Easting values were used as
-        # the stage
-
-        #assert 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 allclose(-elevation, y)  #Meters
-        
-        fid.close()
-        self.delete_mux(files)
-        os.remove(sww_file)
-        
-  
-    def test_urs2sww_origin(self):
-        tide = 1
-        base_name, files = self.create_mux()
-        urs2sww(base_name
-                , origin=(0,0,0)
-                , mean_stage=tide
-                , remove_nc_files=True,
-                      verbose=self.verbose
-                )
-        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)
-
-        #  x and y are absolute
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        geo_reference = Geo_reference(NetCDFObject=fid)
-
-        
-        time = fid.variables['time'][:]
-        #print "time", time
-        assert num.allclose([0.,0.5,1.], time)
-        assert fid.starttime == 0.0
-        #Check that first coordinate is correctly represented       
-        #Work out the UTM coordinates for first point
-        zone, e, n = redfearn(-34.5, 150.66667)       
-       
-        assert num.allclose([x[0],y[0]], [e,n])
-
-        
-        #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)
-        #momentum = velocity*(stage+elevation)
-        # -(-elevation) since elevation is inverted in mux files
-        # = n*(e+tide+n) based on how I'm writing these files
-        answer = n*(e+tide+n)
-        actual = xmomentum[0,0]
-        assert num.allclose(answer, actual)  #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_urs2sww_minmaxlatlong(self):
-        
-        #longitudes = [150.66667, 150.83334, 151., 151.16667]
-        #latitudes = [-34.5, -34.33333, -34.16667, -34]
-
-        tide = 1
-        base_name, files = self.create_mux()
-        urs2sww(base_name,
-                minlat=-34.5,
-                maxlat=-34,
-                minlon= 150.66667,
-                maxlon= 151.16667,
-                mean_stage=tide,
-                remove_nc_files=True,
-                      verbose=self.verbose
-                )
-        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(-34.5, 150.66667) 
-        assert num.allclose([x[0],y[0]], [e,n])
-
-        
-        #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)
-        #momentum = velocity*(stage+elevation)
-        # -(-elevation) since elevation is inverted in mux files
-        # = n*(e+tide+n) based on how I'm writing these files
-        answer = n*(e+tide+n)
-        actual = xmomentum[0,0]
-        assert num.allclose(answer, actual)  #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_urs2sww_minmaxmintmaxt(self):
-        
-        #longitudes = [150.66667, 150.83334, 151., 151.16667]
-        #latitudes = [-34.5, -34.33333, -34.16667, -34]
-
-        tide = 1
-        base_name, files = self.create_mux()
-        
-        urs2sww(base_name,
-                mint=0.25,
-                maxt=0.75,
-                mean_stage=tide,
-                remove_nc_files=True,
-                verbose=self.verbose)
-        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)
-
-        
-        time = fid.variables['time'][:]
-        assert num.allclose(time, [0.0]) # the time is relative
-        assert fid.starttime == 0.5
-        
-        fid.close()
-        self.delete_mux(files)
-        #print "sww_file", sww_file
-        os.remove(sww_file)
-  
-
-#-------------------------------------------------------------
-
-if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_File_Conversion,'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
-
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 7768)
+++ /trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 7769)
@@ -1,5 +1,5 @@
 #!/usr/bin/env python
 
-import unittest, os
+import unittest, os, time
 import os.path
 from math import pi, sqrt
@@ -7,4 +7,7 @@
 
 from anuga.shallow_water import Domain
+
+from Scientific.IO.NetCDF import NetCDFFile
+from anuga.file.sww import extent_sww
 
 from anuga.config import g, epsilon
@@ -14,5 +17,6 @@
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.geospatial_data.geospatial_data import Geospatial_data
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross, \
+                                            rectangular
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
 from anuga.shallow_water.forcing import Inflow, Cross_section
@@ -1465,8 +1469,4 @@
         This run tries it with georeferencing and with elevation = -1
         """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-        from mesh_factory import rectangular
 
         # Create basic mesh (20m x 3m)
@@ -5262,8 +5262,4 @@
         was in March-April 2007
         """
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-        from data_manager import extent_sww
-        from mesh_factory import rectangular
 
         # Create basic mesh
