Changeset 7759

Timestamp:

May 31, 2010, 4:58:56 PM (14 years ago)

Author:

hudson

Message:

Cut out file conversion routines from data_manager.

Location:

trunk/anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (6 diffs)
• data_manager_joaquims_patch.py (deleted)
• test_data_manager.py (modified) (4 diffs)

trunk/anuga_core/source/anuga/shallow_water/data_manager.py

@@ -365 +365 @@
     infile.close()
     outfile.close()
-
-
-
-##
-# @brief  Return block of surface lines for each cross section
-# @param lines Iterble  of text lines to process.
-# @note BROKEN?  UNUSED?
-def _read_hecras_cross_sections(lines):
-    """Return block of surface lines for each cross section
-    Starts with SURFACE LINE,
-    Ends with END CROSS-SECTION
-    """
-
-    points = []
-
-    reading_surface = False
-    for i, line in enumerate(lines):
-        if len(line.strip()) == 0:    #Ignore blanks
-            continue
-
-        if lines[i].strip().startswith('SURFACE LINE'):
-            reading_surface = True
-            continue
-
-        if lines[i].strip().startswith('END') and reading_surface:
-            yield points
-            reading_surface = False
-            points = []
-
-        if reading_surface:
-            fields = line.strip().split(',')
-            easting = float(fields[0])
-            northing = float(fields[1])
-            elevation = float(fields[2])
-            points.append([easting, northing, elevation])
-
-
-##
-# @brief Convert HECRAS (.sdf) file to PTS file.
-# @param basename_in Sterm of input filename.
-# @param basename_out Sterm of output filename.
-# @param verbose True if this function is to be verbose.
-def hecras_cross_sections2pts(basename_in,
-                              basename_out=None,
-                              verbose=False):
-    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
-
-    Example:
-
-# RAS export file created on Mon 15Aug2005 11:42
-# by HEC-RAS Version 3.1.1
-
-BEGIN HEADER:
-  UNITS: METRIC
-  DTM TYPE: TIN
-  DTM: v:\1\cit\perth_topo\river_tin
-  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
-  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
-  MAP PROJECTION: UTM
-  PROJECTION ZONE: 50
-  DATUM: AGD66
-  VERTICAL DATUM:
-  NUMBER OF REACHES:  19
-  NUMBER OF CROSS-SECTIONS:  14206
-END HEADER:
-
-Only the SURFACE LINE data of the following form will be utilised
-  CROSS-SECTION:
-    STREAM ID:Southern-Wungong
-    REACH ID:Southern-Wungong
-    STATION:19040.*
-    CUT LINE:
-      405548.671603161 , 6438142.7594925
-      405734.536092045 , 6438326.10404912
-      405745.130459356 , 6438331.48627354
-      405813.89633823 , 6438368.6272789
-    SURFACE LINE:
-     405548.67,   6438142.76,   35.37
-     405552.24,   6438146.28,   35.41
-     405554.78,   6438148.78,   35.44
-     405555.80,   6438149.79,   35.44
-     405559.37,   6438153.31,   35.45
-     405560.88,   6438154.81,   35.44
-     405562.93,   6438156.83,   35.42
-     405566.50,   6438160.35,   35.38
-     405566.99,   6438160.83,   35.37
-     ...
-   END CROSS-SECTION
-
-    Convert to NetCDF pts format which is
-
-    points:  (Nx2) float array
-    elevation: N float array
-    """
-
-    import os
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    root = basename_in
-
-    # Get ASCII file
-    infile = open(root + '.sdf', 'r')
-
-    if verbose: log.critical('Reading DEM from %s' % (root + '.sdf'))
-
-    lines = infile.readlines()
-    infile.close()
-
-    if verbose: log.critical('Converting to pts format')
-
-    # Scan through the header, picking up stuff we need.
-    i = 0
-    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
-        i += 1
-
-    assert lines[i].strip().upper() == 'BEGIN HEADER:'
-    i += 1
-
-    assert lines[i].strip().upper().startswith('UNITS:')
-    units = lines[i].strip().split()[1]
-    i += 1
-
-    assert lines[i].strip().upper().startswith('DTM TYPE:')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('DTM:')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('STREAM')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('CROSS')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
-    projection = lines[i].strip().split(':')[1]
-    i += 1
-
-    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
-    zone = int(lines[i].strip().split(':')[1])
-    i += 1
-
-    assert lines[i].strip().upper().startswith('DATUM:')
-    datum = lines[i].strip().split(':')[1]
-    i += 1
-
-    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
-    i += 1
-
-    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
-    number_of_cross_sections = int(lines[i].strip().split(':')[1])
-    i += 1
-
-    # Now read all points
-    points = []
-    elevation = []
-    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
-        for k, entry in enumerate(entries):
-            points.append(entry[:2])
-            elevation.append(entry[2])
-
-    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
-          %(j+1, number_of_cross_sections)
-    assert j+1 == number_of_cross_sections, msg
-
-    # Get output file, write PTS data
-    if basename_out == None:
-        ptsname = root + '.pts'
-    else:
-        ptsname = basename_out + '.pts'
-
-    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
-    geo = Geospatial_data(points, {"elevation":elevation},
-                          verbose=verbose, geo_reference=geo_ref)
-    geo.export_points_file(ptsname)
 
 
@@ -729 +551 @@
 
 
-##
-# @brief Convert SWW file to PTS file (at selected points).
-# @param basename_in Stem name of input SWW file.
-# @param basename_out Stem name of output file.
-# @param data_points If given, points where quantity is to be computed.
-# @param quantity Name (or expression) of existing quantity(s) (def: elevation).
-# @param timestep If given, output quantity at that timestep.
-# @param reduction If given, reduce quantity by this factor.
-# @param NODATA_value The NODATA value (default -9999).
-# @param verbose True if this function is to be verbose.
-# @param origin ??
-def sww2pts(basename_in, basename_out=None,
-            data_points=None,
-            quantity=None,
-            timestep=None,
-            reduction=None,
-            NODATA_value=-9999,
-            verbose=False,
-            origin=None):
-    """Read SWW file and convert to interpolated values at selected points
-
-    The parameter 'quantity' must be the name of an existing quantity or
-    an expression involving existing quantities. The default is 'elevation'.
-
-    if timestep (an index) is given, output quantity at that timestep.
-
-    if reduction is given use that to reduce quantity over all timesteps.
-
-    data_points (Nx2 array) give locations of points where quantity is to 
-    be computed.
-    """
-
-    import sys
-    from anuga.geometry.polygon import inside_polygon, outside_polygon
-    from anuga.abstract_2d_finite_volumes.util import \
-             apply_expression_to_dictionary
-    from anuga.geospatial_data.geospatial_data import Geospatial_data
-
-    if quantity is None:
-        quantity = 'elevation'
-
-    if reduction is None:
-        reduction = max
-
-    if basename_out is None:
-        basename_out = basename_in + '_%s' % quantity
-
-    swwfile = basename_in + '.sww'
-    ptsfile = basename_out + '.pts'
-
-    # Read sww file
-    if verbose: log.critical('Reading from %s' % swwfile)
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(swwfile)
-
-    # Get extent and reference
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    volumes = fid.variables['volumes'][:]
-
-    number_of_timesteps = fid.dimensions['number_of_timesteps']
-    number_of_points = fid.dimensions['number_of_points']
-    if origin is None:
-        # Get geo_reference
-        # sww files don't have to have a geo_ref
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference() # Default georef object
-
-        xllcorner = geo_reference.get_xllcorner()
-        yllcorner = geo_reference.get_yllcorner()
-        zone = geo_reference.get_zone()
-    else:
-        zone = origin[0]
-        xllcorner = origin[1]
-        yllcorner = origin[2]
-
-    # FIXME: Refactor using code from file_function.statistics
-    # Something like print swwstats(swwname)
-    if verbose:
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        times = fid.variables['time'][:]
-        log.critical('------------------------------------------------')
-        log.critical('Statistics of SWW file:')
-        log.critical('  Name: %s' % swwfile)
-        log.critical('  Reference:')
-        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
-        log.critical('    Start time: %f' % fid.starttime[0])
-        log.critical('  Extent:')
-        log.critical('    x [m] in [%f, %f], len(x) == %d'
-                     % (num.min(x), num.max(x), len(x.flat)))
-        log.critical('    y [m] in [%f, %f], len(y) == %d'
-                     % (num.min(y), num.max(y), len(y.flat)))
-        log.critical('    t [s] in [%f, %f], len(t) == %d'
-                     % (min(times), max(times), len(times)))
-        log.critical('  Quantities [SI units]:')
-        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
-            q = fid.variables[name][:].flat
-            log.critical('    %s in [%f, %f]' % (name, min(q), max(q)))
-
-    # Get quantity and reduce if applicable
-    if verbose: log.critical('Processing quantity %s' % quantity)
-
-    # Turn NetCDF objects into numeric arrays
-    quantity_dict = {}
-    for name in fid.variables.keys():
-        quantity_dict[name] = fid.variables[name][:]
-
-    # Convert quantity expression to quantities found in sww file
-    q = apply_expression_to_dictionary(quantity, quantity_dict)
-
-    if len(q.shape) == 2:
-        # q has a time component and needs to be reduced along
-        # the temporal dimension
-        if verbose: log.critical('Reducing quantity %s' % quantity)
-
-        q_reduced = num.zeros(number_of_points, num.float)
-        for k in range(number_of_points):
-            q_reduced[k] = reduction(q[:,k])
-        q = q_reduced
-
-    # Post condition: Now q has dimension: number_of_points
-    assert len(q.shape) == 1
-    assert q.shape[0] == number_of_points
-
-    if verbose:
-        log.critical('Processed values for %s are in [%f, %f]'
-                     % (quantity, min(q), max(q)))
-
-    # Create grid and update xll/yll corner and x,y
-    vertex_points = num.concatenate((x[:, num.newaxis], y[:, num.newaxis]), axis=1)
-    assert len(vertex_points.shape) == 2
-
-    # Interpolate
-    from anuga.fit_interpolate.interpolate import Interpolate
-    interp = Interpolate(vertex_points, volumes, verbose=verbose)
-
-    # Interpolate using quantity values
-    if verbose: log.critical('Interpolating')
-    interpolated_values = interp.interpolate(q, data_points).flatten()
-
-    if verbose:
-        log.critical('Interpolated values are in [%f, %f]'
-                     % (num.min(interpolated_values),
-                        num.max(interpolated_values)))
-
-    # Assign NODATA_value to all points outside bounding polygon
-    # (from interpolation mesh)
-    P = interp.mesh.get_boundary_polygon()
-    outside_indices = outside_polygon(data_points, P, closed=True)
-
-    for i in outside_indices:
-        interpolated_values[i] = NODATA_value
-
-    # Store results
-    G = Geospatial_data(data_points=data_points, attributes=interpolated_values)
-
-    G.export_points_file(ptsfile, absolute = True)
-
-    fid.close()
-
-
 
 ##
@@ -933 +591 @@
 # @param very_verbose
 # @return 
-def sww2domain(filename, boundary=None, t=None,
+def load_sww_as_domain(filename, boundary=None, t=None,
                fail_if_NaN=True, NaN_filler=0,
                verbose=False, very_verbose=False):
@@ -1882 +1540 @@
     return d <= max_distance
 
-################################################################################
-# CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
-################################################################################
-
-WAVEHEIGHT_MUX_LABEL = '-z-mux'
-EAST_VELOCITY_LABEL =  '-e-mux'
-NORTH_VELOCITY_LABEL =  '-n-mux'
-
-##
-# @brief Convert URS file(s) (wave prop) to an SWW file.
-# @param basename_in Stem of the input filenames.
-# @param basename_out Path to the output SWW file.
-# @param verbose True if this function is to be verbose.
-# @param mint
-# @param maxt
-# @param mean_stage
-# @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing).
-# @param hole_points_UTM
-# @param zscale
-# @note Also convert latitude and longitude to UTM. All coordinates are
-#       assumed to be given in the GDA94 datum.
-# @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux'
-#       added for relative height, x-velocity and y-velocity respectively.
-def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
-                      mint=None, maxt=None,
-                      mean_stage=0,
-                      origin=None,
-                      hole_points_UTM=None,
-                      zscale=1):
-    """
-    Convert URS C binary format for wave propagation to
-    sww format native to abstract_2d_finite_volumes.
-
-    Specify only basename_in and read files of the form
-    basefilename-z-mux, basefilename-e-mux and
-    basefilename-n-mux containing relative height,
-    x-velocity and y-velocity, respectively.
-
-    Also convert latitude and longitude to UTM. All coordinates are
-    assumed to be given in the GDA94 datum. The latitude and longitude
-    information is assumed ungridded grid.
-
-    min's and max's: If omitted - full extend is used.
-    To include a value min ans max may equal it.
-    Lat and lon are assumed to be in decimal degrees.
-
-    origin is a 3-tuple with geo referenced
-    UTM coordinates (zone, easting, northing)
-    It will be the origin of the sww file. This shouldn't be used,
-    since all of anuga should be able to handle an arbitary origin.
-    The mux point info is NOT relative to this origin.
-
-    URS C binary format has data organised as TIME, LONGITUDE, LATITUDE
-    which means that latitude is the fastest
-    varying dimension (row major order, so to speak)
-
-    In URS C binary the latitudes and longitudes are in assending order.
-
-    Note, interpolations of the resulting sww file will be different
-    from results of urs2sww.  This is due to the interpolation
-    function used, and the different grid structure between urs2sww
-    and this function.
-
-    Interpolating data that has an underlying gridded source can
-    easily end up with different values, depending on the underlying
-    mesh.
-
-    consider these 4 points
-    50  -50
-
-    0     0
-
-    The grid can be
-     -
-    |\|   A
-     -
-     or;
-      -
-     |/|  B
-      -
-
-    If a point is just below the center of the midpoint, it will have a
-    +ve value in grid A and a -ve value in grid B.
-    """
-
-    from anuga.mesh_engine.mesh_engine import NoTrianglesError
-    from anuga.pmesh.mesh import Mesh
-
-    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
-                basename_in + EAST_VELOCITY_LABEL,
-                basename_in + NORTH_VELOCITY_LABEL]
-    quantities = ['HA','UA','VA']
-
-    # instantiate urs_points of the three mux files.
-    mux = {}
-    for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
-
-    # Could check that the depth is the same. (hashing)
-
-    # handle to a mux file to do depth stuff
-    a_mux = mux[quantities[0]]
-
-    # Convert to utm
-    lat = a_mux.lonlatdep[:,1]
-    long = a_mux.lonlatdep[:,0]
-    points_utm, zone = convert_from_latlon_to_utm(latitudes=lat,
-                                                  longitudes=long)
-
-    elevation = a_mux.lonlatdep[:,2] * -1
-
-    # grid (create a mesh from the selected points)
-    # This mesh has a problem.  Triangles are streched over ungridded areas.
-    # If these areas could be described as holes in pmesh, that would be great.
-
-    # I can't just get the user to selection a point in the middle.
-    # A boundary is needed around these points.
-    # But if the zone of points is obvious enough auto-segment should do
-    # a good boundary.
-    mesh = Mesh()
-    mesh.add_vertices(points_utm)
-    mesh.auto_segment(smooth_indents=True, expand_pinch=True)
-
-    # To try and avoid alpha shape 'hugging' too much
-    mesh.auto_segment(mesh.shape.get_alpha() * 1.1)
-    if hole_points_UTM is not None:
-        point = ensure_absolute(hole_points_UTM)
-        mesh.add_hole(point[0], point[1])
-
-    try:
-        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
-    except NoTrianglesError:
-        # This is a bit of a hack, going in and changing the data structure.
-        mesh.holes = []
-        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
-
-    mesh_dic = mesh.Mesh2MeshList()
-
-    #mesh.export_mesh_file(basename_in + '_168.tsh')
-    #import sys; sys.exit()
-    # These are the times of the mux file
-    mux_times = []
-    for i in range(a_mux.time_step_count):
-        mux_times.append(a_mux.time_step * i)
-    (mux_times_start_i, mux_times_fin_i) = mux2sww_time(mux_times, mint, maxt)
-    times = mux_times[mux_times_start_i:mux_times_fin_i]
-
-    if mux_times_start_i == mux_times_fin_i:
-        # Close the mux files
-        for quantity, file in map(None, quantities, files_in):
-            mux[quantity].close()
-        msg = "Due to mint and maxt there's no time info in the boundary SWW."
-        raise Exception, msg
-
-    # If this raise is removed there is currently no downstream errors
-
-    points_utm=ensure_numeric(points_utm)
-    assert num.alltrue(ensure_numeric(mesh_dic['generatedpointlist'])
-                       == ensure_numeric(points_utm))
-
-    volumes = mesh_dic['generatedtrianglelist']
-
-    # Write sww intro and grid stuff.
-    if basename_out is None:
-        swwname = basename_in + '.sww'
-    else:
-        swwname = basename_out + '.sww'
-
-    if verbose: log.critical('Output to %s' % swwname)
-
-    outfile = NetCDFFile(swwname, netcdf_mode_w)
-
-    # For a different way of doing this, check out tsh2sww
-    # work out sww_times and the index range this covers
-    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
-    sww.store_header(outfile, times, len(volumes), len(points_utm),
-                     verbose=verbose, sww_precision=netcdf_float)
-    outfile.mean_stage = mean_stage
-    outfile.zscale = zscale
-
-    sww.store_triangulation(outfile, points_utm, volumes,
-                            zone,  
-                            new_origin=origin,
-                            verbose=verbose)
-    sww.store_static_quantities(outfile, elevation=elevation)
-
-    if verbose: log.critical('Converting quantities')
-
-    # Read in a time slice from each mux file and write it to the SWW file
-    j = 0
-    for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']):
-        if j >= mux_times_start_i and j < mux_times_fin_i:
-            stage = zscale*ha + mean_stage
-            h = stage - elevation
-            xmomentum = ua*h
-            ymomentum = -1 * va * h # -1 since in mux files south is positive.
-            sww.store_quantities(outfile,
-                                 slice_index=j-mux_times_start_i,
-                                 verbose=verbose,
-                                 stage=stage,
-                                 xmomentum=xmomentum,
-                                 ymomentum=ymomentum,
-                                 sww_precision=num.float)
-        j += 1
-
-    if verbose: sww.verbose_quantities(outfile)
-
-    outfile.close()
-
-    # Do some conversions while writing the sww file
 
 
@@ -2214 +1662 @@
 # @param maxt ??
 # @return ??
-def mux2sww_time(mux_times, mint, maxt):
+def read_time_from_mux(mux_times, mint, maxt):
     """
     """
@@ -2233 +1681 @@
 
 
-##
-# @brief Convert a URS (mux2, wave propagation) file to an STS file.
-# @param basename_in String (or list) of source file stems.
-# @param basename_out Stem of output STS file path.
-# @param weights
-# @param verbose True if this function is to be verbose.
-# @param origin Tuple with with geo-ref UTM coords (zone, easting, northing).
-# @param zone
-# @param mean_stage
-# @param zscale
-# @param ordering_filename Path of a file specifying which mux2 gauge points are
-#                          to be stored.
-# @note Also convert latitude and longitude to UTM. All coordinates are
-#       assumed to be given in the GDA94 datum.
-def urs2sts(basename_in, basename_out=None,
-            weights=None,
-            verbose=False,
-            origin=None,
-            zone=None,
-            central_meridian=None,            
-            mean_stage=0.0,
-            zscale=1.0,
-            ordering_filename=None):
-    """Convert URS mux2 format for wave propagation to sts format
-
-    Also convert latitude and longitude to UTM. All coordinates are
-    assumed to be given in the GDA94 datum
-
-    origin is a 3-tuple with geo referenced
-    UTM coordinates (zone, easting, northing)
-
-    inputs:
-
-    basename_in: list of source file prefixes
-
-        These are combined with the extensions:
-        WAVEHEIGHT_MUX2_LABEL = '-z-mux2' for stage
-        EAST_VELOCITY_MUX2_LABEL = '-e-mux2' xmomentum
-        NORTH_VELOCITY_MUX2_LABEL = '-n-mux2' and ymomentum
-
-        to create a 2D list of mux2 file. The rows are associated with each
-        quantity and must have the above extensions
-        the columns are the list of file prefixes.
-
-    ordering: a .txt file name specifying which mux2 gauge points are
-              to be stored. This is indicated by the index of the gauge
-              in the ordering file.
-
-              ordering file format:
-              1st line:    'index,longitude,latitude\n'
-              other lines: index,longitude,latitude
-
-              If ordering is None or ordering file is empty then
-               all points are taken in the order they
-              appear in the mux2 file.
-
-
-    output:
-      basename_out: name of sts file in which mux2 data is stored.
-      
-      
-      
-    NOTE: South is positive in mux files so sign of y-component of velocity is reverted
-    """
-
-    import os
-    from Scientific.IO.NetCDF import NetCDFFile
-    from types import ListType,StringType
-    from operator import __and__
-
-    if not isinstance(basename_in, ListType):
-        if verbose: log.critical('Reading single source')
-        basename_in = [basename_in]
-
-    # This is the value used in the mux file format to indicate NAN data
-    # FIXME (Ole): This should be changed everywhere to IEEE NAN when
-    #              we upgrade to Numpy
-    NODATA = 99
-
-    # Check that basename is a list of strings
-    if not reduce(__and__, map(lambda z:isinstance(z,StringType), basename_in)):
-        msg= 'basename_in must be a string or list of strings'
-        raise Exception, msg
-
-    # Find the number of sources to be used
-    numSrc = len(basename_in)
-
-    # A weight must be specified for each source
-    if weights is None:
-        # Default is equal weighting
-        weights = num.ones(numSrc, num.float) / numSrc
-    else:
-        weights = ensure_numeric(weights)
-        msg = 'When combining multiple sources must specify a weight for ' \
-              'mux2 source file'
-        assert len(weights) == numSrc, msg
-
-    if verbose: log.critical('Weights used in urs2sts: %s' % str(weights))
-
-    # Check output filename
-    if basename_out is None:
-        msg = 'STS filename must be specified as basename_out ' \
-              'in function urs2sts'
-        raise Exception, msg
-
-    if basename_out.endswith('.sts'):
-        stsname = basename_out
-    else:
-        stsname = basename_out + '.sts'
-
-    # Create input filenames from basenames and check their existence
-    files_in = [[], [], []]
-    for files in basename_in:
-        files_in[0].append(files + WAVEHEIGHT_MUX2_LABEL),
-        files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL)
-        files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL)
-
-    quantities = ['HA','UA','VA'] # Quantity names used in the MUX2 format
-    for i in range(len(quantities)):
-        for file_in in files_in[i]:
-            if (os.access(file_in, os.R_OK) == 0):
-                msg = 'File %s does not exist or is not accessible' % file_in
-                raise IOError, msg
-
-    # Establish permutation array
-    if ordering_filename is not None:
-        if verbose is True: log.critical('Reading ordering file %s'
-                                         % ordering_filename)
-
-        # Read ordering file
-        try:
-            fid = open(ordering_filename, 'r')
-            file_header = fid.readline().split(',')
-            ordering_lines = fid.readlines()
-            fid.close()
-        except:
-            msg = 'Cannot open %s' % ordering_filename
-            raise Exception, msg
-
-        reference_header = 'index, longitude, latitude\n'
-        reference_header_split = reference_header.split(',')
-        for i in range(3):
-            if not file_header[i].strip() == reference_header_split[i].strip():
-                msg = 'File must contain header: ' + reference_header
-                raise Exception, msg
-
-        if len(ordering_lines) < 2:
-            msg = 'File must contain at least two points'
-            raise Exception, msg
-
-        permutation = [int(line.split(',')[0]) for line in ordering_lines]
-        permutation = ensure_numeric(permutation)
-    else:
-        permutation = None
-
-    # Read MUX2 files
-    if (verbose): log.critical('reading mux2 file')
-
-    mux={}
-    for i, quantity in enumerate(quantities):
-        # For each quantity read the associated list of source mux2 file with
-        # extention associated with that quantity
-
-        times, latitudes, longitudes, elevation, mux[quantity], starttime \
-            = read_mux2_py(files_in[i], weights, permutation, verbose=verbose)
-
-        # Check that all quantities have consistent time and space information
-        if quantity != quantities[0]:
-            msg = '%s, %s and %s have inconsistent gauge data' \
-                  % (files_in[0], files_in[1], files_in[2])
-            assert num.allclose(times, times_old), msg
-            assert num.allclose(latitudes, latitudes_old), msg
-            assert num.allclose(longitudes, longitudes_old), msg
-            assert num.allclose(elevation, elevation_old), msg
-            assert num.allclose(starttime, starttime_old), msg
-        times_old = times
-        latitudes_old = latitudes
-        longitudes_old = longitudes
-        elevation_old = elevation
-        starttime_old = starttime
-
-        # Self check - can be removed to improve speed
-        #ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]
-        #ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]
-        #
-        #msg = 'Longitudes specified in ordering file do not match those ' \
-        #      'found in mux files. ' \
-        #      'I got %s instead of %s (only beginning shown)' \
-        #      % (str(longitudes[:10]) + '...',
-        #         str(ref_longitudes[:10]) + '...')
-        #assert allclose(longitudes, ref_longitudes), msg
-        #
-        #msg = 'Latitudes specified in ordering file do not match those ' \
-        #      'found in mux files. '
-        #      'I got %s instead of %s (only beginning shown)' \
-        #      % (str(latitudes[:10]) + '...',
-        #         str(ref_latitudes[:10]) + '...')
-        #assert allclose(latitudes, ref_latitudes), msg
-
-    # Store timeseries in STS file
-    msg = 'File is empty and or clipped region not in file region'
-    assert len(latitudes > 0), msg
-
-    number_of_points = latitudes.shape[0]      # Number of stations retrieved
-    number_of_times = times.shape[0]           # Number of timesteps
-    number_of_latitudes = latitudes.shape[0]   # Number latitudes
-    number_of_longitudes = longitudes.shape[0] # Number longitudes
-
-    # The permutation vector of contains original indices
-    # as given in ordering file or None in which case points
-    # are assigned the trivial indices enumerating them from
-    # 0 to number_of_points-1
-    if permutation is None:
-        permutation = num.arange(number_of_points, dtype=num.int)
-
-    # NetCDF file definition
-    outfile = NetCDFFile(stsname, netcdf_mode_w)
-
-    description = 'Converted from URS mux2 files: %s' % basename_in
-
-    # Create new file
-    sts = Write_sts()
-    sts.store_header(outfile,
-                     times+starttime,
-                     number_of_points,
-                     description=description,
-                     verbose=verbose,
-                     sts_precision=netcdf_float)
-
-    # Store
-    from anuga.coordinate_transforms.redfearn import redfearn
-
-    x = num.zeros(number_of_points, num.float)  # Easting
-    y = num.zeros(number_of_points, num.float)  # Northing
-
-    # Check zone boundaries
-    if zone is None:
-        refzone, _, _ = redfearn(latitudes[0], longitudes[0],
-                                 central_meridian=central_meridian)
-    else:
-        refzone = zone
-
-    old_zone = refzone
-
-    for i in range(number_of_points):
-        computed_zone, easting, northing = redfearn(latitudes[i], longitudes[i],
-                                                    zone=zone,
-                                                    central_meridian=central_meridian)
-        x[i] = easting
-        y[i] = northing
-        if computed_zone != refzone:
-            msg = 'All sts gauges need to be in the same zone. \n'
-            msg += 'offending gauge:Zone %d,%.4f, %4f\n' \
-                   % (computed_zone, easting, northing)
-            msg += 'previous gauge:Zone %d,%.4f, %4f' \
-                   % (old_zone, old_easting, old_northing)
-            raise Exception, msg
-        old_zone = computed_zone
-        old_easting = easting
-        old_northing = northing
-
-    if origin is None:
-        origin = Geo_reference(refzone, min(x), min(y))
-    geo_ref = write_NetCDF_georeference(origin, outfile)
-
-    elevation = num.resize(elevation, outfile.variables['elevation'][:].shape)
-    outfile.variables['permutation'][:] = permutation.astype(num.int32) # Opteron 64
-    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
-    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
-    outfile.variables['elevation'][:] = elevation
-
-    stage = outfile.variables['stage']
-    xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum']
-
-    if verbose: log.critical('Converting quantities')
-
-    for j in range(len(times)):
-        for i in range(number_of_points):
-            ha = mux['HA'][i,j]
-            ua = mux['UA'][i,j]
-            va = mux['VA'][i,j]
-            if ha == NODATA:
-                if verbose:
-                    msg = 'Setting nodata value %d to 0 at time = %f, ' \
-                          'point = %d' % (ha, times[j], i)
-                    log.critical(msg)
-                ha = 0.0
-                ua = 0.0
-                va = 0.0
-
-            w = zscale*ha + mean_stage
-            h = w - elevation[i]
-            stage[j,i] = w
-
-            xmomentum[j,i] = ua * h
-            ymomentum[j,i] = -va * h # South is positive in mux files
-
-
-    outfile.close()
-
-
-##
-# @brief Create a list of points defining a boundary from an STS file.
-# @param stsname Stem of path to the STS file to read.
-# @return A list of boundary points.
-def create_sts_boundary(stsname):
-    """Create a list of points defining a boundary from an STS file.
-
-    Create boundary segments from .sts file. Points can be stored in
-    arbitrary order within the .sts file. The order in which the .sts points
-    make up the boundary are given in order.txt file
-
-    FIXME:
-    Point coordinates are stored in relative eastings and northings.
-    But boundary is produced in absolute coordinates
-    """
-
-    try:
-        fid = NetCDFFile(stsname + '.sts', netcdf_mode_r)
-    except:
-        msg = 'Cannot open %s' % stsname + '.sts'
-        raise msg
-
-    xllcorner = fid.xllcorner[0]
-    yllcorner = fid.yllcorner[0]
-
-    #Points stored in sts file are normalised to [xllcorner,yllcorner] but
-    #we cannot assume that boundary polygon will be. At least the
-    #additional points specified by the user after this function is called
-    x = fid.variables['x'][:] + xllcorner
-    y = fid.variables['y'][:] + yllcorner
-
-    x = num.reshape(x, (len(x),1))
-    y = num.reshape(y, (len(y),1))
-    sts_points = num.concatenate((x,y), axis=1)
-
-    return sts_points.tolist()
-
-
-##
-# @brief Obsolete?
-# @param outfile
-# @param has
-# @param uas
-# @param vas
-# @param elevation
-# @param mean_stage
-# @param zscale
-# @param verbose
-def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation,
-                                   mean_stage=0, zscale=1,
-                                   verbose=False):
-    #Time stepping
-    stage = outfile.variables['stage']
-    xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum']
-
-    n = len(has)
-    j = 0
-    for ha, ua, va in map(None, has, uas, vas):
-        if verbose and j % ((n+10)/10) == 0: log.critical('  Doing %d of %d'
-                                                          % (j, n))
-        w = zscale*ha + mean_stage
-        stage[j] = w
-        h = w - elevation
-        xmomentum[j] = ua * h
-        ymomentum[j] = -1 * va * h  # -1 since in mux files south is positive.
-        j += 1
-
-
-##
-# @brief Convert a set of URS files to a text file.
-# @param basename_in Stem path to the 3 URS files.
-# @param location_index ??
-def urs2txt(basename_in, location_index=None):
-    """
-    Not finished or tested
-    """
-
-    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
-                basename_in + EAST_VELOCITY_LABEL,
-                basename_in + NORTH_VELOCITY_LABEL]
-    quantities = ['HA','UA','VA']
-
-    d = ","
-
-    # instantiate urs_points of the three mux files.
-    mux = {}
-    for quantity, file in map(None, quantities, files_in):
-        mux[quantity] = Urs_points(file)
-
-    # Could check that the depth is the same. (hashing)
-
-    # handle to a mux file to do depth stuff
-    a_mux = mux[quantities[0]]
-
-    # Convert to utm
-    latitudes = a_mux.lonlatdep[:,1]
-    longitudes = a_mux.lonlatdep[:,0]
-    points_utm, zone = \
-        convert_from_latlon_to_utm(latitudes=latitudes, longitudes=longitudes)
-    depths = a_mux.lonlatdep[:,2]
-
-    # open the output text file, start writing.
-    fid = open(basename_in + '.txt', 'w')
-
-    fid.write("zone: " + str(zone) + "\n")
-
-    if location_index is not None:
-        #Title
-        li = location_index
-        fid.write('location_index' + d + 'lat' + d + 'long' + d +
-                  'Easting' + d + 'Northing' + '\n')
-        fid.write(str(li) + d + str(latitudes[li]) + d +
-                  str(longitudes[li]) + d + str(points_utm[li][0]) + d +
-                  str(points_utm[li][01]) + '\n')
-
-    # the non-time dependent stuff
-    #Title
-    fid.write('location_index' + d + 'lat' + d + 'long' + d +
-              'Easting' + d + 'Northing' + d + 'depth m' + '\n')
-    i = 0
-    for depth, point_utm, lat, long in map(None, depths, points_utm,
-                                           latitudes, longitudes):
-
-        fid.write(str(i) + d + str(lat) + d + str(long) + d +
-                  str(point_utm[0]) + d + str(point_utm[01]) + d +
-                  str(depth) + '\n')
-        i += 1
-
-    #Time dependent
-    if location_index is not None:
-        time_step = a_mux.time_step
-        i = 0
-        #Title
-        fid.write('time' + d + 'HA depth m' + d + 'UA momentum East x m/sec' +
-                  d + 'VA momentum North y m/sec' + '\n')
-        for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
-            fid.write(str(i*time_step) + d + str(HA[location_index]) + d +
-                      str(UA[location_index]) + d +
-                      str(VA[location_index]) + '\n')
-            i += 1
-
-
-
-
-##
-# @brief A class to write STS files.
-class Write_sts:
-    sts_quantities = ['stage','xmomentum','ymomentum']
-    RANGE = '_range'
-    EXTREMA = ':extrema'
-
-    ##
-    # @brief Instantiate this STS writer class.
-    def __init__(self):
-        pass
-
-    ##
-    # @brief Write a header to the underlying data file.
-    # @param outfile Handle to open file to write.
-    # @param times A list of the time slice times *or* a start time.
-    # @param number_of_points The number of URS gauge sites.
-    # @param description Description string to write into the STS file.
-    # @param sts_precision Format of data to write (netcdf constant ONLY).
-    # @param verbose True if this function is to be verbose.
-    # @note If 'times' is a list, the info will be made relative.
-    def store_header(self,
-                     outfile,
-                     times,
-                     number_of_points,
-                     description='Converted from URS mux2 format',
-                     sts_precision=netcdf_float32,
-                     verbose=False):
-        """
-        outfile - the name of the file that will be written
-        times - A list of the time slice times OR a start time
-        Note, if a list is given the info will be made relative.
-        number_of_points - the number of urs gauges sites
-        """
-
-        outfile.institution = 'Geoscience Australia'
-        outfile.description = description
-
-        try:
-            revision_number = get_revision_number()
-        except:
-            revision_number = None
-
-        # Allow None to be stored as a string
-        outfile.revision_number = str(revision_number)
-
-        # Start time in seconds since the epoch (midnight 1/1/1970)
-        # This is being used to seperate one number from a list.
-        # what it is actually doing is sorting lists from numeric arrays.
-        if isinstance(times, (list, num.ndarray)):
-            number_of_times = len(times)
-            times = ensure_numeric(times)
-            if number_of_times == 0:
-                starttime = 0
-            else:
-                starttime = times[0]
-                times = times - starttime  #Store relative times
-        else:
-            number_of_times = 0
-            starttime = times
-
-        outfile.starttime = starttime
-
-        # Dimension definitions
-        outfile.createDimension('number_of_points', number_of_points)
-        outfile.createDimension('number_of_timesteps', number_of_times)
-        outfile.createDimension('numbers_in_range', 2)
-
-        # Variable definitions
-        outfile.createVariable('permutation', netcdf_int, ('number_of_points',))
-        outfile.createVariable('x', sts_precision, ('number_of_points',))
-        outfile.createVariable('y', sts_precision, ('number_of_points',))
-        outfile.createVariable('elevation',sts_precision, ('number_of_points',))
-
-        q = 'elevation'
-        outfile.createVariable(q + Write_sts.RANGE, sts_precision,
-                               ('numbers_in_range',))
-
-        # Initialise ranges with small and large sentinels.
-        # If this was in pure Python we could have used None sensibly
-        outfile.variables[q + Write_sts.RANGE][0] = max_float  # Min
-        outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
-
-        # Doing sts_precision instead of Float gives cast errors.
-        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
-
-        for q in Write_sts.sts_quantities:
-            outfile.createVariable(q, sts_precision, ('number_of_timesteps',
-                                                      'number_of_points'))
-            outfile.createVariable(q + Write_sts.RANGE, sts_precision,
-                                   ('numbers_in_range',))
-            # Initialise ranges with small and large sentinels.
-            # If this was in pure Python we could have used None sensibly
-            outfile.variables[q + Write_sts.RANGE][0] = max_float  # Min
-            outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
-
-        if isinstance(times, (list, num.ndarray)):
-            outfile.variables['time'][:] = times    #Store time relative
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('Statistics:')
-            log.critical('    t in [%f, %f], len(t) == %d'
-                         % (num.min(times), num.max(times), len(times.flat)))
-
-    ##
-    # @brief 
-    # @param outfile 
-    # @param points_utm 
-    # @param elevation 
-    # @param zone 
-    # @param new_origin 
-    # @param points_georeference 
-    # @param verbose True if this function is to be verbose.
-    def store_points(self,
-                     outfile,
-                     points_utm,
-                     elevation, zone=None, new_origin=None,
-                     points_georeference=None, verbose=False):
-
-        """
-        points_utm - currently a list or array of the points in UTM.
-        points_georeference - the georeference of the points_utm
-
-        How about passing new_origin and current_origin.
-        If you get both, do a convertion from the old to the new.
-
-        If you only get new_origin, the points are absolute,
-        convert to relative
-
-        if you only get the current_origin the points are relative, store
-        as relative.
-
-        if you get no georefs create a new georef based on the minimums of
-        points_utm.  (Another option would be to default to absolute)
-
-        Yes, and this is done in another part of the code.
-        Probably geospatial.
-
-        If you don't supply either geo_refs, then supply a zone. If not
-        the default zone will be used.
-
-        precondition:
-             header has been called.
-        """
-
-        number_of_points = len(points_utm)
-        points_utm = num.array(points_utm)
-
-        # given the two geo_refs and the points, do the stuff
-        # described in the method header
-        points_georeference = ensure_geo_reference(points_georeference)
-        new_origin = ensure_geo_reference(new_origin)
-
-        if new_origin is None and points_georeference is not None:
-            points = points_utm
-            geo_ref = points_georeference
-        else:
-            if new_origin is None:
-                new_origin = Geo_reference(zone, min(points_utm[:,0]),
-                                                 min(points_utm[:,1]))
-            points = new_origin.change_points_geo_ref(points_utm,
-                                                      points_georeference)
-            geo_ref = new_origin
-
-        # At this stage I need a georef and points
-        # the points are relative to the georef
-        geo_ref.write_NetCDF(outfile)
-
-        x = points[:,0]
-        y = points[:,1]
-        z = outfile.variables['elevation'][:]
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('More Statistics:')
-            log.critical('  Extent (/lon):')
-            log.critical('    x in [%f, %f], len(lat) == %d'
-                         % (min(x), max(x), len(x)))
-            log.critical('    y in [%f, %f], len(lon) == %d'
-                         % (min(y), max(y), len(y)))
-            log.critical('    z in [%f, %f], len(z) == %d'
-                         % (min(elevation), max(elevation), len(elevation)))
-            log.critical('geo_ref: %s' % str(geo_ref))
-            log.critical('------------------------------------------------')
-
-        z = resize(bath_grid,outfile.variables['elevation'][:].shape)
-        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
-        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
-        #outfile.variables['z'][:] = elevation
-        outfile.variables['elevation'][:] = elevation  #FIXME HACK4
-
-        # This updates the _range values
-        q = 'elevation'
-        outfile.variables[q + Write_sts.RANGE][0] = min(elevation)
-        outfile.variables[q + Write_sts.RANGE][1] = max(elevation)
-
-    ##
-    # @brief Store quantity data in the underlying file.
-    # @param outfile
-    # @param sts_precision
-    # @param slice_index
-    # @param time
-    # @param verboseTrue if this function is to be verbose.
-    # @param **quant Extra keyword args.
-    def store_quantities(self, outfile, sts_precision=num.float32,
-                         slice_index=None, time=None,
-                         verbose=False, **quant):
-        """Write the quantity info.
-
-        **quant is extra keyword arguments passed in. These must be
-          the sts quantities, currently; stage.
-
-        if the time array is already been built, use the slice_index
-        to specify the index.
-
-        Otherwise, use time to increase the time dimension
-
-        Maybe make this general, but the viewer assumes these quantities,
-        so maybe we don't want it general - unless the viewer is general
-
-        precondition:
-            triangulation and header have been called.
-        """
-
-        if time is not None:
-            file_time = outfile.variables['time']
-            slice_index = len(file_time)
-            file_time[slice_index] = time
-
-        # Write the conserved quantities from Domain.
-        # Typically stage, xmomentum, ymomentum
-        # other quantities will be ignored, silently.
-        # Also write the ranges: stage_range
-        for q in Write_sts.sts_quantities:
-            if not quant.has_key(q):
-                msg = 'STS file can not write quantity %s' % q
-                raise NewQuantity, msg
-            else:
-                q_values = quant[q]
-                outfile.variables[q][slice_index] = \
-                                q_values.astype(sts_precision)
-
-                # This updates the _range values
-                q_range = outfile.variables[q + Write_sts.RANGE][:]
-                q_values_min = num.min(q_values)
-                if q_values_min < q_range[0]:
-                    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
-                q_values_max = num.max(q_values)
-                if q_values_max > q_range[1]:
-                    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
-
-
-##
-# @brief 
-class Urs_points:
-    """
-    Read the info in URS mux files.
-
-    for the quantities here's a correlation between the file names and
-    what they mean;
-    z-mux is height above sea level, m
-    e-mux is velocity is Eastern direction, m/s
-    n-mux is velocity is Northern direction, m/s
-    """
-
-    ##
-    # @brief Initialize this instance of Urs_points.
-    # @param urs_file Path to the underlying data file.
-    def __init__(self, urs_file):
-        self.iterated = False
-        columns = 3                         # long, lat , depth
-        mux_file = open(urs_file, 'rb')
-
-        # Number of points/stations
-        (self.points_num,) = unpack('i', mux_file.read(4))
-
-        # nt, int - Number of time steps
-        (self.time_step_count,) = unpack('i', mux_file.read(4))
-        #dt, float - time step, seconds
-        (self.time_step,) = unpack('f', mux_file.read(4))
-        msg = "Bad data in the urs file."
-        if self.points_num < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-        if self.time_step_count < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-        if self.time_step < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-
-        # The depth is in meters, and it is the distance from the ocean
-        # to the sea bottom.
-        lonlatdep = p_array.array('f')
-        lonlatdep.read(mux_file, columns * self.points_num)
-        lonlatdep = num.array(lonlatdep, dtype=num.float)
-        lonlatdep = num.reshape(lonlatdep, (self.points_num, columns))
-        self.lonlatdep = lonlatdep
-
-        self.mux_file = mux_file
-        # check this array
-
-    ##
-    # @brief Allow iteration over quantity data wrt time.
-    def __iter__(self):
-        """
-        iterate over quantity data which is with respect to time.
-
-        Note: You can only iterate once over an object
-
-        returns quantity infomation for each time slice
-        """
-
-        msg =  "You can only interate once over a urs file."
-        assert not self.iterated, msg
-
-        self.iter_time_step = 0
-        self.iterated = True
-
-        return self
-
-    ##
-    # @brief 
-    def next(self):
-        if self.time_step_count == self.iter_time_step:
-            self.close()
-            raise StopIteration
-
-        #Read in a time slice from mux file
-        hz_p_array = p_array.array('f')
-        hz_p_array.read(self.mux_file, self.points_num)
-        hz_p = num.array(hz_p_array, dtype=num.float)
-        self.iter_time_step += 1
-
-        return hz_p
-
-    ##
-    # @brief Close the mux file.
-    def close(self):
-        self.mux_file.close()
-
-################################################################################
-# END URS UNGRIDDED 2 SWW
-################################################################################
+
 
 

trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -18 +18 @@
 from anuga.shallow_water.data_manager import *
 from anuga.shallow_water.sww_file import SWW_file
-from anuga.file_conversion.file_conversion import tsh2sww, \
-                        asc_csiro2sww, pmesh_to_domain_instance
 from anuga.coordinate_transforms.geo_reference import Geo_reference                        
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
@@ -738 +736 @@
 
 
-    def test_hecras_cross_sections2pts(self):
-        """Test conversion from HECRAS cross sections in ascii format
-        to native NetCDF pts format
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Write test asc file
-        root = 'hecrastest'
-
-        filename = root+'.sdf'
-        fid = open(filename, 'w')
-        fid.write("""
-# RAS export file created on Mon 15Aug2005 11:42
-# by HEC-RAS Version 3.1.1
-
-BEGIN HEADER:
-  UNITS: METRIC
-  DTM TYPE: TIN
-  DTM: v:\1\cit\perth_topo\river_tin
-  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
-  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
-  MAP PROJECTION: UTM
-  PROJECTION ZONE: 50
-  DATUM: AGD66
-  VERTICAL DATUM:
-  NUMBER OF REACHES:  19
-  NUMBER OF CROSS-SECTIONS:  2
-END HEADER:
-
-
-BEGIN CROSS-SECTIONS:
-
-  CROSS-SECTION:
-    STREAM ID:Southern-Wungong
-    REACH ID:Southern-Wungong
-    STATION:21410
-    CUT LINE:
-      407546.08 , 6437277.542
-      407329.32 , 6437489.482
-      407283.11 , 6437541.232
-    SURFACE LINE:
-     407546.08,   6437277.54,   52.14
-     407538.88,   6437284.58,   51.07
-     407531.68,   6437291.62,   50.56
-     407524.48,   6437298.66,   49.58
-     407517.28,   6437305.70,   49.09
-     407510.08,   6437312.74,   48.76
-  END:
-
-  CROSS-SECTION:
-    STREAM ID:Swan River
-    REACH ID:Swan Mouth
-    STATION:840.*
-    CUT LINE:
-      381178.0855 , 6452559.0685
-      380485.4755 , 6453169.272
-    SURFACE LINE:
-     381178.09,   6452559.07,   4.17
-     381169.49,   6452566.64,   4.26
-     381157.78,   6452576.96,   4.34
-     381155.97,   6452578.56,   4.35
-     381143.72,   6452589.35,   4.43
-     381136.69,   6452595.54,   4.58
-     381114.74,   6452614.88,   4.41
-     381075.53,   6452649.43,   4.17
-     381071.47,   6452653.00,   3.99
-     381063.46,   6452660.06,   3.67
-     381054.41,   6452668.03,   3.67
-  END:
-END CROSS-SECTIONS:
-""")
-
-        fid.close()
-
-
-        #Convert to NetCDF pts
-        hecras_cross_sections2pts(root)
-
-        #Check contents
-        #Get NetCDF
-        fid = NetCDFFile(root+'.pts', netcdf_mode_r)
-
-        # Get the variables
-        #print fid.variables.keys()
-        points = fid.variables['points']
-        elevation = fid.variables['elevation']
-
-        #Check values
-        ref_points = [[407546.08, 6437277.54],
-                      [407538.88, 6437284.58],
-                      [407531.68, 6437291.62],
-                      [407524.48, 6437298.66],
-                      [407517.28, 6437305.70],
-                      [407510.08, 6437312.74]]
-
-        ref_points += [[381178.09, 6452559.07],
-                       [381169.49, 6452566.64],
-                       [381157.78, 6452576.96],
-                       [381155.97, 6452578.56],
-                       [381143.72, 6452589.35],
-                       [381136.69, 6452595.54],
-                       [381114.74, 6452614.88],
-                       [381075.53, 6452649.43],
-                       [381071.47, 6452653.00],
-                       [381063.46, 6452660.06],
-                       [381054.41, 6452668.03]]
-
-
-        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
-        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
-
-        #print points[:]
-        #print ref_points
-        assert num.allclose(points, ref_points)
-
-        #print attributes[:]
-        #print ref_elevation
-        assert num.allclose(elevation, ref_elevation)
-
-        #Cleanup
-        fid.close()
-
-
-        os.remove(root + '.sdf')
-        os.remove(root + '.pts')
-
 
     def test_export_grid(self):
@@ -1359 +1229 @@
         os.remove(ascfile)
         os.remove(swwfile2)
-
-
-
-    def test_sww2pts_centroids(self):
-        """Test that sww information can be converted correctly to pts data at specified coordinates
-        - in this case, the centroids.
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-        # Used for points that lie outside mesh
-        NODATA_value = 1758323
-
-        # Setup
-        self.domain.set_name('datatest')
-
-        ptsfile = self.domain.get_name() + '_elevation.pts'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.smooth = True #self.set_store_vertices_uniquely(False)
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        # Check contents in NetCDF
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        elevation = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-        volumes = fid.variables['volumes'][:]
-
-
-        # Invoke interpolation for vertex points       
-        points = num.concatenate( (x[:,num.newaxis],y[:,num.newaxis]), axis=1 )
-        points = num.ascontiguousarray(points)
-        sww2pts(self.domain.get_name(),
-                quantity = 'elevation',
-                data_points = points,
-                NODATA_value = NODATA_value,
-                verbose = self.verbose)
-        ref_point_values = elevation
-        point_values = Geospatial_data(ptsfile).get_attributes()
-        #print 'P', point_values
-        #print 'Ref', ref_point_values        
-        assert num.allclose(point_values, ref_point_values)        
-
-
-
-        # Invoke interpolation for centroids
-        points = self.domain.get_centroid_coordinates()
-        #print points
-        sww2pts(self.domain.get_name(),
-                quantity = 'elevation',
-                data_points = points,
-                NODATA_value = NODATA_value,
-                verbose = self.verbose)
-        ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids
-
-        
-        point_values = Geospatial_data(ptsfile).get_attributes()
-        #print 'P', point_values
-        #print 'Ref', ref_point_values        
-        assert num.allclose(point_values, ref_point_values)        
-
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(sww.filename)
-        os.remove(ptsfile)
 
 
@@ -2027 +1811 @@
 
         os.remove(root + '.dem')
-        os.remove(root + '_100.dem')
-
-
-    def write_mux(self, lat_long_points, time_step_count, time_step,
-                  depth=None, ha=None, ua=None, va=None):
-        """
-        This will write 3 non-gridded mux files, for testing.
-        If no quantities are passed in,
-        na and va quantities will be the Easting values.
-        Depth and ua will be the Northing value.
-        
-        The mux file format has south as positive so
-        this function will swap the sign for va.  
-        """
-
-        #print "lat_long_points", lat_long_points
-        #print "time_step_count",time_step_count
-        #print "time_step",
-
-        
-        points_num = len(lat_long_points)
-        lonlatdeps = []
-        quantities = ['HA','UA','VA']
-        
-        mux_names = [WAVEHEIGHT_MUX_LABEL,
-                     EAST_VELOCITY_LABEL,
-                     NORTH_VELOCITY_LABEL]
-        quantities_init = [[],[],[]]
-        # urs binary is latitude fastest
-        for point in lat_long_points:
-            lat = point[0]
-            lon = point[1]
-            _ , e, n = redfearn(lat, lon)
-            if depth is None:
-                this_depth = n
-            else:
-                this_depth = depth
-            if ha is None:
-                this_ha = e
-            else:
-                this_ha = ha
-            if ua is None:
-                this_ua = n
-            else:
-                this_ua = ua
-            if va is None:
-                this_va = e   
-            else:
-                this_va = va         
-            lonlatdeps.append([lon, lat, this_depth])
-            quantities_init[0].append(this_ha) # HA
-            quantities_init[1].append(this_ua) # UA
-            quantities_init[2].append(this_va) # VA 
-                
-        file_handle, base_name = tempfile.mkstemp("")
-        os.close(file_handle)
-        os.remove(base_name)
-
-        files = []        
-        for i, q in enumerate(quantities): 
-            quantities_init[i] = ensure_numeric(quantities_init[i])
-            #print "HA_init", HA_init
-            q_time = num.zeros((time_step_count, points_num), num.float)
-            for time in range(time_step_count):
-                q_time[time,:] = quantities_init[i] #* time * 4
-            
-            #Write C files
-            columns = 3 # long, lat , depth
-            file = base_name + mux_names[i]
-            #print "base_name file",file 
-            f = open(file, 'wb')
-            files.append(file)
-            f.write(pack('i',points_num))
-            f.write(pack('i',time_step_count))
-            f.write(pack('f',time_step))
-
-            #write lat/long info
-            for lonlatdep in lonlatdeps:
-                for float in lonlatdep:
-                    f.write(pack('f',float))
-                    
-            # Write quantity info
-            for time in  range(time_step_count):
-                for point_i in range(points_num):
-                    f.write(pack('f',q_time[time,point_i]))
-                    #print " mux_names[i]", mux_names[i] 
-                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
-            f.close()
-        return base_name, files
-        
-    
-    def delete_mux(self, files):
-        for file in files:
-            os.remove(file)
-        
-    def write_mux2(self, lat_long_points, time_step_count, time_step,
-                   first_tstep, last_tstep,
-                   depth=None, ha=None, ua=None, va=None):
-        """
-        This will write 3 non-gridded mux files, for testing.
-        If no quantities are passed in,
-        na and va quantities will be the Easting values.
-        Depth and ua will be the Northing value.
-        """
-        #print "lat_long_points", lat_long_points
-        #print "time_step_count",time_step_count
-        #print "time_step",
-
-        #irrelevant header information
-        ig=ilon=ilat=0
-        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
-
-        points_num = len(lat_long_points)
-        latlondeps = []
-        quantities = ['HA','UA','VA']
-
-        mux_names = [WAVEHEIGHT_MUX2_LABEL,
-                     EAST_VELOCITY_MUX2_LABEL,
-                     NORTH_VELOCITY_MUX2_LABEL]
-
-        msg='first_tstep and last_step arrays must have same length as number of points'
-        assert len(first_tstep)==points_num,msg
-        assert len(last_tstep)==points_num,msg
-
-        if depth is not None:
-            depth=ensure_numeric(depth)
-            assert len(depth)==points_num
-        if ha is not None:
-            ha=ensure_numeric(ha)
-            assert ha.shape==(points_num,time_step_count)
-        if ua is not None:
-            ua=ensure_numeric(ua)
-            assert ua.shape==(points_num,time_step_count)
-        if va is not None:
-            va=ensure_numeric(va)
-            assert va.shape==(points_num,time_step_count)
-
-        quantities_init = [[],[],[]]
-        # urs binary is latitude fastest
-        for i,point in enumerate(lat_long_points):
-            lat = point[0]
-            lon = point[1]
-            _ , e, n = redfearn(lat, lon)
-            if depth is None:
-                this_depth = n
-            else:
-                this_depth = depth[i]
-            latlondeps.append([lat, lon, this_depth])
-
-            if ha is None:
-                this_ha = e
-                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
-            else:
-                quantities_init[0].append(ha[i])
-            if ua is None:
-                this_ua = n
-                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
-            else:
-                quantities_init[1].append(ua[i])
-            if va is None:
-                this_va = e
-                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
-            else:
-                quantities_init[2].append(-va[i]) # South is negative in MUX
-
-        file_handle, base_name = tempfile.mkstemp("write_mux2")
-        os.close(file_handle)
-        os.remove(base_name)
-
-        files = []        
-        for i, q in enumerate(quantities):
-            q_time = num.zeros((time_step_count, points_num), num.float)
-            quantities_init[i] = ensure_numeric(quantities_init[i])
-            for time in range(time_step_count):
-                #print i, q, time, quantities_init[i][:,time]
-                q_time[time,:] = quantities_init[i][:,time]
-                #print i, q, time, q_time[time, :]                
-
-            #Write C files
-            columns = 3 # long, lat , depth
-            file = base_name + mux_names[i]
-            
-            #print 'base_name file', file 
-            f = open(file, 'wb')
-            files.append(file)
-
-            f.write(pack('i',points_num))
-            #write mux 2 header
-            for latlondep in latlondeps:
-                f.write(pack('f',latlondep[0]))
-                f.write(pack('f',latlondep[1]))
-                f.write(pack('f',mcolat))
-                f.write(pack('f',mcolon))
-                f.write(pack('i',ig))
-                f.write(pack('i',ilon))
-                f.write(pack('i',ilat))
-                f.write(pack('f',latlondep[2]))
-                f.write(pack('f',centerlat))
-                f.write(pack('f',centerlon))
-                f.write(pack('f',offset))
-                f.write(pack('f',az))
-                f.write(pack('f',baz))
-                f.write(pack('f',time_step))
-                f.write(pack('i',time_step_count))
-                for j in range(4): # identifier
-                    f.write(pack('f',id))    
-
-            #first_tstep=1
-            #last_tstep=time_step_count
-            for i,latlondep in enumerate(latlondeps):
-                f.write(pack('i',first_tstep[i]))
-            for i,latlondep in enumerate(latlondeps):
-                f.write(pack('i',last_tstep[i]))
-
-            # Find when first station starts recording
-            min_tstep = min(first_tstep)
-            # Find when all stations have stopped recording
-            max_tstep = max(last_tstep)
-
-            #for time in  range(time_step_count):
-            for time in range(min_tstep-1,max_tstep):
-                f.write(pack('f',time*time_step))                
-                for point_i in range(points_num):
-                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
-                        #print 'writing', time, point_i, q_time[time, point_i]
-                        f.write(pack('f', q_time[time, point_i]))
-            f.close()
-
-        return base_name, files
-
-    def test_urs2sts_read_mux2_pyI(self):
-        """test_urs2sts_read_mux2_pyI(self):
-        Constant stage,momentum at each gauge
-        """
-        tide = 1
-        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)
-        last_tstep=time_step_count*num.ones(n,num.int)
-        depth=20*num.ones(n,num.float)
-        ha=2*num.ones((n,time_step_count),num.float)
-        ua=5*num.ones((n,time_step_count),num.float)
-        va=-10*num.ones((n,time_step_count),num.float)
-        #-ve added to take into account mux file format where south is positive.
-        base_name, files = self.write_mux2(lat_long_points,
-                                      time_step_count, time_step,
-                                      first_tstep, last_tstep,
-                                      depth=depth,
-                                      ha=ha,
-                                      ua=ua,
-                                      va=va)
-
-        weights=num.ones(1, num.float)
-        #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage, starttime = read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
-        msg='ha and ua have different gauge meta data'
-        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity, starttime_va=read_mux2_py([files[2]],weights)
-        msg='ha and va have different gauge meta data'
-        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
-
-        self.delete_mux(files)
-
-        msg='time array has incorrect length'
-        assert times.shape[0]==time_step_count,msg
-        
-        msg = 'time array is incorrect'
-        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
-        assert num.allclose(times,time_step*num.arange(time_step_count)), msg
-        
-        msg='Incorrect gauge positions returned'
-        for i,point in enumerate(lat_long_points):
-            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
-
-        msg='Incorrect gauge depths returned'
-        assert num.allclose(elevation,-depth),msg
-        msg='incorrect gauge height time series returned'
-        assert num.allclose(stage,ha)
-        msg='incorrect gauge ua time series returned'
-        assert num.allclose(xvelocity,ua)
-        msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity, -va)
-
-    def test_urs2sts_read_mux2_pyII(self):
-        """Spatially varing stage
-        """
-        tide = 1
-        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)
-        last_tstep=(time_step_count)*num.ones(n,num.int)
-        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)+1
-        ha[1]=time_step_count-num.arange(1,time_step_count+1)
-        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)
-        #-ve added to take into account mux file format where south is positive.
-        base_name, files = self.write_mux2(lat_long_points,
-                                      time_step_count, time_step,
-                                      first_tstep, last_tstep,
-                                      depth=depth,
-                                      ha=ha,
-                                      ua=ua,
-                                      va=va)
-
-        weights=num.ones(1, num.float)
-        #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
-        msg='ha and ua have different gauge meta data'
-        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
-        msg='ha and va have different gauge meta data'
-        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
-
-
-        self.delete_mux(files)
-
-        msg='time array has incorrect length'
-        #assert times.shape[0]==time_step_count,msg
-        msg = 'time array is incorrect'
-        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
-        msg='Incorrect gauge positions returned'
-        for i,point in enumerate(lat_long_points):
-            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
-
-        msg='Incorrect gauge depths returned'
-        assert num.allclose(elevation, -depth),msg
-        msg='incorrect gauge height time series returned'
-        assert num.allclose(stage, ha)
-        msg='incorrect gauge ua time series returned'
-        assert num.allclose(xvelocity, ua)
-        msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity, -va) # South is positive in MUX
-
-    def test_urs2sts_read_mux2_pyIII(self):
-        """Varying start and finish times
-        """
-        tide = 1
-        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
-
-        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)
-        #-ve added to take into account mux file format where south is positive.
-        base_name, files = self.write_mux2(lat_long_points,
-                                      time_step_count, time_step,
-                                      first_tstep, last_tstep,
-                                      depth=depth,
-                                      ha=ha,
-                                      ua=ua,
-                                      va=va)
-
-        weights=num.ones(1, num.float)
-        #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage, starttime=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity, starttime_ua=read_mux2_py([files[1]],weights)
-        msg='ha and ua have different gauge meta data'
-        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
-        msg='ha and va have different gauge meta data'
-        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
-
-        self.delete_mux(files)
-
-        msg='time array has incorrect length'
-        #assert times.shape[0]==time_step_count,msg
-        msg = 'time array is incorrect'
-        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
-        msg='Incorrect gauge positions returned'
-        for i,point in enumerate(lat_long_points):
-            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
-
-
-        # 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;
-        msg='Incorrect gauge depths returned'
-        assert num.allclose(elevation,-depth),msg
-        msg='incorrect gauge height time series returned'
-        assert num.allclose(stage,ha)
-        msg='incorrect gauge ua time series returned'
-        assert num.allclose(xvelocity,ua)
-        msg='incorrect gauge va time series returned'
-        assert num.allclose(yvelocity, -va) # South is positive in mux
-        
-
-        
-    def test_read_mux_platform_problem1(self):
-        """test_read_mux_platform_problem1
-        
-        This is to test a situation where read_mux returned 
-        wrong values Win32
-
-        This test passes on Windows but test_read_mux_platform_problem2
-        does not
-        """
-        
-        from urs_ext import read_mux2 
-        
-        verbose = False
-                
-        tide = 1.5
-        time_step_count = 10
-        time_step = 0.2
-        times_ref = num.arange(0, time_step_count*time_step, time_step)
-
-        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
-        n = len(lat_long_points)
-        
-        # 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        
-        
-        # 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
-            
-        # 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                                  
-        
-        # 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 ordering file
-        permutation = ensure_numeric([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 files back and check values
-        weights = ensure_numeric([1.0])
-
-        # For each quantity read the associated list of source mux2 file with 
-        # extention associated with that quantity
-        file_params=-1*num.ones(3,num.float) #[nsta,dt,nt]
-        OFFSET = 5
-
-        for j, file in enumerate(filesI):
-            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
-
-            number_of_selected_stations = data.shape[0]
-
-            # Index where data ends and parameters begin
-            parameters_index = data.shape[1]-OFFSET          
-          
-            for i in range(number_of_selected_stations):
-                if j == 0: assert num.allclose(data[i][:parameters_index], ha0[permutation[i], :])
-                if j == 1: assert num.allclose(data[i][:parameters_index], ua0[permutation[i], :])
-                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
-        
-        self.delete_mux(filesI)
-
-
-        
-        
-    def test_read_mux_platform_problem2(self):
-        """test_read_mux_platform_problem2
-        
-        This is to test a situation where read_mux returned 
-        wrong values Win32
-
-        This test does not pass on Windows but test_read_mux_platform_problem1
-        does
-        """
-        
-        from urs_ext import read_mux2 
-        
-        from anuga.config import single_precision as epsilon        
-        
-        verbose = False
-                
-        tide = 1.5
-        time_step_count = 10
-        time_step = 0.2
-        
-        times_ref = num.arange(0, time_step_count*time_step, time_step)
-        
-        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
-                           (-21.,115.), (-22., 117.)]
-        n = len(lat_long_points)
-        
-        # 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        
-        
-        # 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
-            
-        # 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)
-        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 'Second station to be written to MUX'
-        #print 'ha', ha1[0,:]
-        #print 'ua', ua1[0,:]
-        #print 'va', va1[0,:]
-        
-        # 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)
-
-
-
-
-        # Read mux file back and verify it's correcness
-
-        ####################################################
-        # FIXME (Ole): This is where the test should
-        # verify that the MUX files are correct.
-
-        #JJ: It appears as though
-        #that certain quantities are not being stored with enough precision
-        #inn muxfile or more likely that they are being cast into a
-        #lower precision when read in using read_mux2 Time step and q_time
-        # are equal but only to approx 1e-7
-        ####################################################
-
-        #define information as it should be stored in mus2 files
-        points_num=len(lat_long_points)
-        depth=gauge_depth
-        ha=ha1
-        ua=ua1
-        va=va1
-        
-        quantities = ['HA','UA','VA']
-        mux_names = [WAVEHEIGHT_MUX2_LABEL,
-                     EAST_VELOCITY_MUX2_LABEL,
-                     NORTH_VELOCITY_MUX2_LABEL]
-        quantities_init = [[],[],[]]
-        latlondeps = []
-        #irrelevant header information
-        ig=ilon=ilat=0
-        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
-        # urs binary is latitude fastest
-        for i,point in enumerate(lat_long_points):
-            lat = point[0]
-            lon = point[1]
-            _ , e, n = redfearn(lat, lon)
-            if depth is None:
-                this_depth = n
-            else:
-                this_depth = depth[i]
-            latlondeps.append([lat, lon, this_depth])
-
-            if ha is None:
-                this_ha = e
-                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
-            else:
-                quantities_init[0].append(ha[i])
-            if ua is None:
-                this_ua = n
-                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
-            else:
-                quantities_init[1].append(ua[i])
-            if va is None:
-                this_va = e
-                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
-            else:
-                quantities_init[2].append(va[i])
-
-        for i, q in enumerate(quantities):
-            #print
-            #print i, q
-            
-            q_time = num.zeros((time_step_count, points_num), num.float)
-            quantities_init[i] = ensure_numeric(quantities_init[i])
-            for time in range(time_step_count):
-                #print i, q, time, quantities_init[i][:,time]
-                q_time[time,:] = quantities_init[i][:,time]
-                #print i, q, time, q_time[time, :]
-
-            
-            filename = base_nameII + mux_names[i]
-            f = open(filename, 'rb')
-            if self.verbose: print 'Reading' + filename
-            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
-            #write mux 2 header
-            for latlondep in latlondeps:
-                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
-                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
-                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
-                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
-                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
-                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
-                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
-                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
-                assert abs(az-unpack('f',f.read(4))[0])<epsilon
-                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
-                
-                x = unpack('f', f.read(4))[0]
-                #print time_step
-                #print x
-                assert abs(time_step-x)<epsilon
-                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
-                for j in range(4): # identifier
-                    assert abs(id-unpack('i',f.read(4))[0])<epsilon 
-
-            #first_tstep=1
-            #last_tstep=time_step_count
-            for i,latlondep in enumerate(latlondeps):
-                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
-            for i,latlondep in enumerate(latlondeps):
-                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
-
-            # Find when first station starts recording
-            min_tstep = min(first_tstep)
-            # Find when all stations have stopped recording
-            max_tstep = max(last_tstep)
-
-            #for time in  range(time_step_count):
-            for time in range(min_tstep-1,max_tstep):
-                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
-                for point_i in range(points_num):
-                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
-                        x = unpack('f',f.read(4))[0]
-                        #print time, x, q_time[time, point_i]
-                        if q == 'VA': x = -x # South is positive in MUX
-                        assert abs(q_time[time, point_i]-x)<epsilon
-
-            f.close()
-                                               
-        # Create ordering file
-        permutation = ensure_numeric([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 files back and check values
-        weights = ensure_numeric([1.0])
-
-        # For each quantity read the associated list of source mux2 file with 
-        # extention associated with that quantity
-        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
-        OFFSET = 5
-
-        for j, file in enumerate(filesII):
-            # Read stage, u, v enumerated as j
-            #print 'Reading', j, file
-            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
-
-            #print 'Data received by Python'
-            #print data[1][8]
-            number_of_selected_stations = data.shape[0]
-
-            # Index where data ends and parameters begin
-            parameters_index = data.shape[1]-OFFSET          
-                 
-            quantity=num.zeros((number_of_selected_stations, parameters_index), num.float)
-            
-            
-            for i in range(number_of_selected_stations):
-        
-                #print i, parameters_index
-                #print quantity[i][:]
-                if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
-                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
-                if j == 2:
-                    # FIXME (Ole): This is where the output is wrong on Win32
-                    
-                    #print
-                    #print j, i
-                    #print 'Input'
-                    #print 'u', ua1[permutation[i], 8]       
-                    #print 'v', va1[permutation[i], 8]
-                
-                    #print 'Output'
-                    #print 'v ', data[i][:parameters_index][8]  
-
-                    # South is positive in MUX
-                    #print "data[i][:parameters_index]", data[i][:parameters_index]
-                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
-                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
-        
-        self.delete_mux(filesII)
-           
-    def test_read_mux_platform_problem3(self):
-        
-        # This is to test a situation where read_mux returned 
-        # wrong values Win32
-
-        
-        from urs_ext import read_mux2 
-        
-        from anuga.config import single_precision as epsilon        
-        
-        verbose = False
-                
-        tide = 1.5
-        time_step_count = 10
-        time_step = 0.02
-
-        '''
-        Win results
-        time_step = 0.2000001
-        This is OK        
-        '''
-        
-        '''
-        Win results
-        time_step = 0.20000001
-
-        ======================================================================
-ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
-----------------------------------------------------------------------
-Traceback (most recent call last):
-  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
-    ha1[0]=num.sin(times_ref)
-ValueError: matrices are not aligned for copy
-
-        '''
-
-        '''
-        Win results
-        time_step = 0.200000001
-        FAIL
-         assert num.allclose(data[i][:parameters_index],
-         -va1[permutation[i], :])
-        '''
-        times_ref = num.arange(0, time_step_count*time_step, time_step)
-        #print "times_ref", times_ref
-        
-        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
-                           (-21.,115.), (-22., 117.)]
-        stations = len(lat_long_points)
-        
-        # Create different timeseries starting and ending at different times 
-        first_tstep=num.ones(stations, 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(stations, 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        
-        
-        # Create varying elevation data (positive values for seafloor)
-        gauge_depth=20*num.ones(stations, num.float)
-        for i in range(stations):
-            gauge_depth[i] += i**2
-            
-        # Create data to be written to second mux file        
-        ha1=num.ones((stations,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((stations,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((stations,time_step_count),num.float)
-        va1[0]=2*num.cos(times_ref-0.87)        
-        va1[1]=3*num.ones(time_step_count)
-        va1[3]=2*num.sin(times_ref-0.71)        
-        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
-        # Ensure data used to write mux file to be zero when gauges are
-        # not recording
-        for i in range(stations):
-             # 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 'Second station to be written to MUX'
-        #print 'ha', ha1[0,:]
-        #print 'ua', ua1[0,:]
-        #print 'va', va1[0,:]
-        
-        # 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)
-        #print "filesII", filesII
-
-
-
-
-        # Read mux file back and verify it's correcness
-
-        ####################################################
-        # FIXME (Ole): This is where the test should
-        # verify that the MUX files are correct.
-
-        #JJ: It appears as though
-        #that certain quantities are not being stored with enough precision
-        #inn muxfile or more likely that they are being cast into a
-        #lower precision when read in using read_mux2 Time step and q_time
-        # are equal but only to approx 1e-7
-        ####################################################
-
-        #define information as it should be stored in mus2 files
-        points_num=len(lat_long_points)
-        depth=gauge_depth
-        ha=ha1
-        ua=ua1
-        va=va1
-        
-        quantities = ['HA','UA','VA']
-        mux_names = [WAVEHEIGHT_MUX2_LABEL,
-                     EAST_VELOCITY_MUX2_LABEL,
-                     NORTH_VELOCITY_MUX2_LABEL]
-        quantities_init = [[],[],[]]
-        latlondeps = []
-        #irrelevant header information
-        ig=ilon=ilat=0
-        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
-        # urs binary is latitude fastest
-        for i,point in enumerate(lat_long_points):
-            lat = point[0]
-            lon = point[1]
-            _ , e, n = redfearn(lat, lon)
-            if depth is None:
-                this_depth = n
-            else:
-                this_depth = depth[i]
-            latlondeps.append([lat, lon, this_depth])
-
-            if ha is None:
-                this_ha = e
-                quantities_init[0].append(num.ones(time_step_count,
-                                                   num.float)*this_ha) # HA
-            else:
-                quantities_init[0].append(ha[i])
-            if ua is None:
-                this_ua = n
-                quantities_init[1].append(num.ones(time_step_count,
-                                                   num.float)*this_ua) # UA
-            else:
-                quantities_init[1].append(ua[i])
-            if va is None:
-                this_va = e
-                quantities_init[2].append(num.ones(time_step_count,
-                                                   num.float)*this_va) #
-            else:
-                quantities_init[2].append(va[i])
-
-        for i, q in enumerate(quantities):
-            #print
-            #print i, q
-            
-            q_time = num.zeros((time_step_count, points_num), num.float)
-            quantities_init[i] = ensure_numeric(quantities_init[i])
-            for time in range(time_step_count):
-                #print i, q, time, quantities_init[i][:,time]
-                q_time[time,:] = quantities_init[i][:,time]
-                #print i, q, time, q_time[time, :]
-
-            
-            filename = base_nameII + mux_names[i]
-            f = open(filename, 'rb')
-            if self.verbose: print 'Reading' + filename
-            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
-            #write mux 2 header
-            for latlondep in latlondeps:
-                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
-                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
-                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
-                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
-                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
-                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
-                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
-                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
-                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
-                assert abs(az-unpack('f',f.read(4))[0])<epsilon
-                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
-                
-                x = unpack('f', f.read(4))[0]
-                #print time_step
-                #print x
-                assert abs(time_step-x)<epsilon
-                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
-                for j in range(4): # identifier
-                    assert abs(id-unpack('i',f.read(4))[0])<epsilon 
-
-            #first_tstep=1
-            #last_tstep=time_step_count
-            for i,latlondep in enumerate(latlondeps):
-                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
-            for i,latlondep in enumerate(latlondeps):
-                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
-
-            # Find when first station starts recording
-            min_tstep = min(first_tstep)
-            # Find when all stations have stopped recording
-            max_tstep = max(last_tstep)
-
-            #for time in  range(time_step_count):
-            for time in range(min_tstep-1,max_tstep):
-                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
-                for point_i in range(points_num):
-                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
-                        x = unpack('f',f.read(4))[0]
-                        #print time, x, q_time[time, point_i]
-                        if q == 'VA': x = -x # South is positive in MUX
-                        #print q+" q_time[%d, %d] = %f" %(time, point_i, 
-                                                      #q_time[time, point_i])
-                        assert abs(q_time[time, point_i]-x)<epsilon
-
-            f.close()
-                            
-        permutation = ensure_numeric([4,0,2])
-                   
-        # Create ordering file
-#         _, 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 files back and check values
-        weights = ensure_numeric([1.0])
-
-        # For each quantity read the associated list of source mux2 file with 
-        # extention associated with that quantity
-        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
-        OFFSET = 5
-
-        for j, file in enumerate(filesII):
-            # Read stage, u, v enumerated as j
-            #print 'Reading', j, file
-            #print "file", file
-            data = read_mux2(1, [file], weights, file_params,
-                             permutation, verbose)
-            #print str(j) + "data", data
-
-            #print 'Data received by Python'
-            #print data[1][8]
-            number_of_selected_stations = data.shape[0]
-            #print "number_of_selected_stations", number_of_selected_stations
-            #print "stations", stations
-
-            # Index where data ends and parameters begin
-            parameters_index = data.shape[1]-OFFSET          
-                 
-            for i in range(number_of_selected_stations):
-        
-                #print i, parameters_index
-                if j == 0:
-                    assert num.allclose(data[i][:parameters_index],
-                                        ha1[permutation[i], :])
-                    
-                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
-                if j == 2:
-                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
-        
-        self.delete_mux(filesII)            
+        os.remove(root + '_100.dem')      
         
     def test_urs2sts0(self):