Index: /anuga_core/source/anuga/file_conversion/__init__.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/__init__.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/__init__.py	(revision 7742)
@@ -0,0 +1,26 @@
+"""Make directory available as a Python package
+"""
+
+# Add path of package to PYTHONPATH to allow C-extensions to be loaded
+import sys
+sys.path += __path__
+
+# Make selected classes available directly
+from boundaries import Reflective_boundary,\
+     Transmissive_momentum_set_stage_boundary,\
+     Dirichlet_discharge_boundary,\
+     Field_boundary,\
+     Transmissive_stage_zero_momentum_boundary,\
+     Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
+
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary, Dirichlet_boundary, \
+            Time_boundary, File_boundary, AWI_boundary
+
+from shallow_water_domain import Domain
+
+
+#from shallow_water_balanced_domain import Swb_domain
+
+
+ 
Index: /anuga_core/source/anuga/file_conversion/dem2pts.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/dem2pts.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/dem2pts.py	(revision 7742)
@@ -0,0 +1,250 @@
+# external modules
+import numpy as num
+
+# ANUGA modules
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+                            
+##
+# @brief Convert DEM data  to PTS data.
+# @param basename_in Stem of input filename.
+# @param basename_out Stem of output filename.
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
+# @param use_cache 
+# @param verbose 
+# @return 
+def dem2pts(basename_in, basename_out=None,
+            easting_min=None, easting_max=None,
+            northing_min=None, northing_max=None,
+            use_cache=False, verbose=False,):
+    """Read Digitial Elevation model from the following NetCDF format (.dem)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Convert to NetCDF pts format which is
+
+    points:  (Nx2) float array
+    elevation: N float array
+    """
+
+    kwargs = {'basename_out': basename_out,
+              'easting_min': easting_min,
+              'easting_max': easting_max,
+              'northing_min': northing_min,
+              'northing_max': northing_max,
+              'verbose': verbose}
+
+    if use_cache is True:
+        from caching import cache
+        result = cache(_dem2pts, basename_in, kwargs,
+                       dependencies = [basename_in + '.dem'],
+                       verbose = verbose)
+
+    else:
+        result = apply(_dem2pts, [basename_in], kwargs)
+
+    return result
+
+
+##
+# @brief 
+# @param basename_in 
+# @param basename_out 
+# @param verbose 
+# @param easting_min 
+# @param easting_max 
+# @param northing_min 
+# @param northing_max 
+def _dem2pts(basename_in, basename_out=None, verbose=False,
+            easting_min=None, easting_max=None,
+            northing_min=None, northing_max=None):
+    """Read Digitial Elevation model from the following NetCDF format (.dem)
+
+    Internal function. See public function dem2pts for details.
+    """
+
+    # FIXME: Can this be written feasibly using write_pts?
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    root = basename_in
+
+    # Get NetCDF
+    infile = NetCDFFile(root + '.dem', netcdf_mode_r) 
+
+    if verbose: log.critical('Reading DEM from %s' % (root + '.dem'))
+
+    ncols = infile.ncols[0]
+    nrows = infile.nrows[0]
+    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
+    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
+    cellsize = infile.cellsize[0]
+    NODATA_value = infile.NODATA_value[0]
+    dem_elevation = infile.variables['elevation']
+
+    zone = infile.zone[0]
+    false_easting = infile.false_easting[0]
+    false_northing = infile.false_northing[0]
+
+    # Text strings
+    projection = infile.projection
+    datum = infile.datum
+    units = infile.units
+
+    # Get output file
+    if basename_out == None:
+        ptsname = root + '.pts'
+    else:
+        ptsname = basename_out + '.pts'
+
+    if verbose: log.critical('Store to NetCDF file %s' % ptsname)
+
+    # NetCDF file definition
+    outfile = NetCDFFile(ptsname, netcdf_mode_w)
+
+    # Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF pts format for compact and portable ' \
+                          'storage of spatial point data'
+
+    # Assign default values
+    if easting_min is None: easting_min = xllcorner
+    if easting_max is None: easting_max = xllcorner + ncols*cellsize
+    if northing_min is None: northing_min = yllcorner
+    if northing_max is None: northing_max = yllcorner + nrows*cellsize
+
+    # Compute offsets to update georeferencing
+    easting_offset = xllcorner - easting_min
+    northing_offset = yllcorner - northing_min
+
+    # Georeferencing
+    outfile.zone = zone
+    outfile.xllcorner = easting_min # Easting of lower left corner
+    outfile.yllcorner = northing_min # Northing of lower left corner
+    outfile.false_easting = false_easting
+    outfile.false_northing = false_northing
+
+    outfile.projection = projection
+    outfile.datum = datum
+    outfile.units = units
+
+    # Grid info (FIXME: probably not going to be used, but heck)
+    outfile.ncols = ncols
+    outfile.nrows = nrows
+
+    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
+    totalnopoints = nrows*ncols
+
+    # Calculating number of NODATA_values for each row in clipped region
+    # FIXME: use array operations to do faster
+    nn = 0
+    k = 0
+    i1_0 = 0
+    j1_0 = 0
+    thisj = 0
+    thisi = 0
+    for i in range(nrows):
+        y = (nrows-i-1)*cellsize + yllcorner
+        for j in range(ncols):
+            x = j*cellsize + xllcorner
+            if easting_min <= x <= easting_max \
+               and northing_min <= y <= northing_max:
+                thisj = j
+                thisi = i
+                if dem_elevation_r[i,j] == NODATA_value:
+                    nn += 1
+
+                if k == 0:
+                    i1_0 = i
+                    j1_0 = j
+
+                k += 1
+
+    index1 = j1_0
+    index2 = thisj
+
+    # Dimension definitions
+    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
+    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
+
+    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
+    nopoints = clippednopoints-nn
+
+    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
+
+    if verbose:
+        log.critical('There are %d values in the elevation' % totalnopoints)
+        log.critical('There are %d values in the clipped elevation'
+                     % clippednopoints)
+        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
+
+    outfile.createDimension('number_of_points', nopoints)
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+
+    # Variable definitions
+    outfile.createVariable('points', netcdf_float, ('number_of_points',
+                                                    'number_of_dimensions'))
+    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
+
+    # Get handles to the variables
+    points = outfile.variables['points']
+    elevation = outfile.variables['elevation']
+
+    lenv = index2-index1+1
+
+    # Store data
+    global_index = 0
+    # for i in range(nrows):
+    for i in range(i1_0, thisi+1, 1):
+        if verbose and i % ((nrows+10)/10) == 0:
+            log.critical('Processing row %d of %d' % (i, nrows))
+
+        lower_index = global_index
+
+        v = dem_elevation_r[i,index1:index2+1]
+        no_NODATA = num.sum(v == NODATA_value)
+        if no_NODATA > 0:
+            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
+        else:
+            newcols = lenv              # ncols_in_bounding_box
+
+        telev = num.zeros(newcols, num.float)
+        tpoints = num.zeros((newcols, 2), num.float)
+
+        local_index = 0
+
+        y = (nrows-i-1)*cellsize + yllcorner
+        #for j in range(ncols):
+        for j in range(j1_0,index2+1,1):
+            x = j*cellsize + xllcorner
+            if easting_min <= x <= easting_max \
+               and northing_min <= y <= northing_max \
+               and dem_elevation_r[i,j] != NODATA_value:
+                tpoints[local_index, :] = [x-easting_min, y-northing_min]
+                telev[local_index] = dem_elevation_r[i, j]
+                global_index += 1
+                local_index += 1
+
+        upper_index = global_index
+
+        if upper_index == lower_index + newcols:
+            points[lower_index:upper_index, :] = tpoints
+            elevation[lower_index:upper_index] = telev
+
+    assert global_index == nopoints, 'index not equal to number of points'
+
+    infile.close()
+    outfile.close()
+
Index: /anuga_core/source/anuga/file_conversion/ferret2sww.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/ferret2sww.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/ferret2sww.py	(revision 7742)
@@ -0,0 +1,433 @@
+"""
+    Convert a ferret file to an SWW file.
+"""
+# external modules
+import numpy as num
+
+
+# ANUGA modules
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_float
+from anuga.shallow_water.sww_file import Write_sww  
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     write_NetCDF_georeference                          
+import anuga.utilities.log as log
+
+#local modules
+from file_conversion import get_min_max_indexes                            
+
+
+def ferret2sww(basename_in, basename_out=None,
+               verbose=False,
+               minlat=None, maxlat=None,
+               minlon=None, maxlon=None,
+               mint=None, maxt=None, mean_stage=0,
+               origin=None, zscale=1,
+               fail_on_NaN=True,
+               NaN_filler=0,
+               elevation=None,
+               inverted_bathymetry=True
+               ): #FIXME: Bathymetry should be obtained
+                                  #from MOST somehow.
+                                  #Alternatively from elsewhere
+                                  #or, as a last resort,
+                                  #specified here.
+                                  #The value of -100 will work
+                                  #for the Wollongong tsunami
+                                  #scenario but is very hacky
+    """Convert MOST and 'Ferret' NetCDF format for wave propagation to
+    sww format native to abstract_2d_finite_volumes.
+
+    Specify only basename_in and read files of the form
+    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc 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.
+
+    min's and max's: If omitted - full extend is used.
+    To include a value min may equal it, while max must exceed it.
+    Lat and lon are assuemd to be in decimal degrees
+
+    origin is a 3-tuple with geo referenced
+    UTM coordinates (zone, easting, northing)
+
+    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
+    which means that longitude is the fastest
+    varying dimension (row major order, so to speak)
+
+    ferret2sww uses grid points as vertices in a triangular grid
+    counting vertices from lower left corner upwards, then right
+    """
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    _assert_lat_long(minlat, maxlat, minlon, maxlon)
+
+    # Get NetCDF data
+    if verbose: log.critical('Reading files %s_*.nc' % basename_in)
+
+    # Wave amplitude (cm)
+    file_h = NetCDFFile(basename_in + '_ha.nc', netcdf_mode_r) 
+    
+    # Velocity (x) (cm/s)
+    file_u = NetCDFFile(basename_in + '_ua.nc', netcdf_mode_r)
+     
+    # Velocity (y) (cm/s)
+    file_v = NetCDFFile(basename_in + '_va.nc', netcdf_mode_r)
+    
+    # Elevation (z) (m)
+    file_e = NetCDFFile(basename_in + '_e.nc', netcdf_mode_r)  
+
+    if basename_out is None:
+        swwname = basename_in + '.sww'
+    else:
+        swwname = basename_out + '.sww'
+
+    # Get dimensions of file_h
+    for dimension in file_h.dimensions.keys():
+        if dimension[:3] == 'LON':
+            dim_h_longitude = dimension
+        if dimension[:3] == 'LAT':
+            dim_h_latitude = dimension
+        if dimension[:4] == 'TIME':
+            dim_h_time = dimension
+
+    times = file_h.variables[dim_h_time]
+    latitudes = file_h.variables[dim_h_latitude]
+    longitudes = file_h.variables[dim_h_longitude]
+
+    kmin, kmax, lmin, lmax = get_min_max_indexes(latitudes[:],
+                                                  longitudes[:],
+                                                  minlat, maxlat,
+                                                  minlon, maxlon)
+    # get dimensions for file_e
+    for dimension in file_e.dimensions.keys():
+        if dimension[:3] == 'LON':
+            dim_e_longitude = dimension
+        if dimension[:3] == 'LAT':
+            dim_e_latitude = dimension
+
+    # get dimensions for file_u
+    for dimension in file_u.dimensions.keys():
+        if dimension[:3] == 'LON':
+            dim_u_longitude = dimension
+        if dimension[:3] == 'LAT':
+            dim_u_latitude = dimension
+
+    # get dimensions for file_v
+    for dimension in file_v.dimensions.keys():
+        if dimension[:3] == 'LON':
+            dim_v_longitude = dimension
+        if dimension[:3] == 'LAT':
+            dim_v_latitude = dimension
+
+    # Precision used by most for lat/lon is 4 or 5 decimals
+    e_lat = num.around(file_e.variables[dim_e_latitude][:], 5)
+    e_lon = num.around(file_e.variables[dim_e_longitude][:], 5)
+
+    # Check that files are compatible
+    assert num.allclose(latitudes, file_u.variables[dim_u_latitude])
+    assert num.allclose(latitudes, file_v.variables[dim_v_latitude])
+    assert num.allclose(latitudes, e_lat)
+
+    assert num.allclose(longitudes, file_u.variables[dim_u_longitude])
+    assert num.allclose(longitudes, file_v.variables[dim_v_longitude])
+    assert num.allclose(longitudes, e_lon)
+
+    if mint is None:
+        jmin = 0
+        mint = times[0]
+    else:
+        jmin = num.searchsorted(times, mint)
+        
+        # numpy.int32 didn't work in slicing of amplitude below
+        jmin = int(jmin)
+
+    if maxt is None:
+        jmax = len(times)
+        maxt = times[-1]
+    else:
+        jmax = num.searchsorted(times, maxt)
+        
+        # numpy.int32 didn't work in slicing of amplitude below
+        jmax = int(jmax)        
+
+    kmin, kmax, lmin, lmax = get_min_max_indexes(latitudes[:],
+                                                  longitudes[:],
+                                                  minlat, maxlat,
+                                                  minlon, maxlon)
+
+
+    times = times[jmin:jmax]
+    latitudes = latitudes[kmin:kmax]
+    longitudes = longitudes[lmin:lmax]
+
+    if verbose: log.critical('cropping')
+
+    zname = 'ELEVATION'
+
+    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
+    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
+    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
+    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
+
+    # Get missing values
+    nan_ha = file_h.variables['HA'].missing_value[0]
+    nan_ua = file_u.variables['UA'].missing_value[0]
+    nan_va = file_v.variables['VA'].missing_value[0]
+    if hasattr(file_e.variables[zname],'missing_value'):
+        nan_e  = file_e.variables[zname].missing_value[0]
+    else:
+        nan_e = None
+
+    # Cleanup
+    missing = (amplitudes == nan_ha)
+    if num.sometrue (missing):
+        if fail_on_NaN:
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_ha.nc'
+            raise DataMissingValuesError, msg
+        else:
+            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
+
+    missing = (uspeed == nan_ua)
+    if num.sometrue (missing):
+        if fail_on_NaN:
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_ua.nc'
+            raise DataMissingValuesError, msg
+        else:
+            uspeed = uspeed*(missing==0) + missing*NaN_filler
+
+    missing = (vspeed == nan_va)
+    if num.sometrue (missing):
+        if fail_on_NaN:
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_va.nc'
+            raise DataMissingValuesError, msg
+        else:
+            vspeed = vspeed*(missing==0) + missing*NaN_filler
+
+    missing = (elevations == nan_e)
+    if num.sometrue (missing):
+        if fail_on_NaN:
+            msg = 'NetCDFFile %s contains missing values' \
+                  % basename_in + '_e.nc'
+            raise DataMissingValuesError, msg
+        else:
+            elevations = elevations*(missing==0) + missing*NaN_filler
+
+    number_of_times = times.shape[0]
+    number_of_latitudes = latitudes.shape[0]
+    number_of_longitudes = longitudes.shape[0]
+
+    assert amplitudes.shape[0] == number_of_times
+    assert amplitudes.shape[1] == number_of_latitudes
+    assert amplitudes.shape[2] == number_of_longitudes
+
+    if verbose:
+        _show_stats((latitudes, longitudes), times, amplitudes, \
+                    (uspeed, vspeed), elevations)
+
+    # print number_of_latitudes, number_of_longitudes
+    number_of_points = number_of_latitudes * number_of_longitudes
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
+
+    file_h.close()
+    file_u.close()
+    file_v.close()
+    file_e.close()
+
+    # NetCDF file definition
+    outfile = NetCDFFile(swwname, netcdf_mode_w)
+
+    description = 'Converted from Ferret files: %s, %s, %s, %s' \
+                  % (basename_in + '_ha.nc',
+                     basename_in + '_ua.nc',
+                     basename_in + '_va.nc',
+                     basename_in + '_e.nc')
+
+    # Create new file
+    starttime = times[0]
+
+    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
+    sww.store_header(outfile, times, number_of_volumes,
+                     number_of_points, description=description,
+                     verbose=verbose, sww_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
+
+    if verbose:
+        log.critical('Making triangular grid')
+
+    # Check zone boundaries
+    refzone, _, _ = redfearn(latitudes[0], longitudes[0])
+
+    vertices = {}
+    i = 0
+    for k, lat in enumerate(latitudes):       # Y direction
+        for l, lon in enumerate(longitudes):  # X direction
+            vertices[l, k] = i
+
+            _, easting, northing = redfearn(lat, lon)
+
+            #msg = 'Zone boundary crossed at longitude =', lon
+            #assert zone == refzone, msg
+            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
+            x[i] = easting
+            y[i] = northing
+            i += 1
+
+    #Construct 2 triangles per 'rectangular' element
+    volumes = []
+    for l in range(number_of_longitudes-1):    # X direction
+        for k in range(number_of_latitudes-1): # Y direction
+            v1 = vertices[l, k+1]
+            v2 = vertices[l, k]
+            v3 = vertices[l+1, k+1]
+            v4 = vertices[l+1, k]
+
+            volumes.append([v1, v2, v3]) #Upper element
+            volumes.append([v4, v3, v2]) #Lower element
+
+    volumes = num.array(volumes, num.int)      #array default#
+
+    if origin is None:
+        origin = Geo_reference(refzone, min(x), min(y))
+    geo_ref = write_NetCDF_georeference(origin, outfile)
+
+    if elevation is not None:
+        z = elevation
+    else:
+        if inverted_bathymetry:
+            z = -1 * elevations
+        else:
+            z = elevations
+    #FIXME: z should be obtained from MOST and passed in here
+
+    #FIXME use the Write_sww instance(sww) to write this info
+    z = num.resize(z, outfile.variables['elevation'][:].shape)
+    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
+    #outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
+    outfile.variables['elevation'][:] = z
+    outfile.variables['volumes'][:] = volumes.astype(num.int32) #For Opteron 64
+
+    #Time stepping
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum = outfile.variables['ymomentum']
+
+    if verbose:
+        log.critical('Converting quantities')
+
+    n = len(times)
+    for j in range(n):
+        if verbose and j % ((n+10)/10) == 0:
+            log.critical('  Doing %d of %d' % (j, n))
+
+        i = 0
+        for k in range(number_of_latitudes):      # Y direction
+            for l in range(number_of_longitudes): # X direction
+                w = zscale * amplitudes[j, k, l] / 100 + mean_stage
+                stage[j, i] = w
+                h = w - z[i]
+                xmomentum[j, i] = uspeed[j, k, l]/100*h
+                ymomentum[j, i] = vspeed[j, k, l]/100*h
+                i += 1
+
+    #outfile.close()
+
+    #FIXME: Refactor using code from file_function.statistics
+    #Something like print swwstats(swwname)
+    if verbose:
+        time_info = times, starttime, mint, maxt
+        _show_sww_stats(outfile, swwname, geo_ref, time_info)
+
+    outfile.close()
+
+
+def _show_stats(latlong, times, amplitudes, speeds, elevations):
+    """ Print the statistics nicely to the log file """
+    
+    latitudes, longitudes = latlong
+    uspeed, vspeed = speeds
+    
+    log.critical('------------------------------------------------')
+    log.critical('Statistics:')
+    log.critical('  Extent (lat/lon):')
+    log.critical('    lat in [%f, %f], len(lat) == %d'
+                 % (num.min(latitudes), num.max(latitudes),
+                    len(latitudes.flat)))
+    log.critical('    lon in [%f, %f], len(lon) == %d'
+                 % (num.min(longitudes), num.max(longitudes),
+                    len(longitudes.flat)))
+    log.critical('    t in [%f, %f], len(t) == %d'
+                 % (num.min(times), num.max(times), len(times.flat)))
+
+    name = 'Amplitudes (ha) [cm]'
+    log.critical('  %s in [%f, %f]'
+                 % (name, num.min(amplitudes), num.max(amplitudes)))
+
+    name = 'Speeds (ua) [cm/s]'
+    log.critical('  %s in [%f, %f]'
+                 % (name, num.min(uspeed), num.max(uspeed)))
+
+    name = 'Speeds (va) [cm/s]'
+    log.critical('  %s in [%f, %f]'
+                 % (name, num.min(vspeed), num.max(vspeed)))
+
+    name = 'Elevations (e) [m]'
+    log.critical('  %s in [%f, %f]'
+                 % (name, num.min(elevations), num.max(elevations)))
+                 
+             
+
+def _show_sww_stats(outfile, swwname, geo_ref, time_info):         
+    """ Log SWW output stats. """
+    times, starttime, mint, maxt = time_info
+    x = outfile.variables['x'][:]
+    y = outfile.variables['y'][:]
+    log.critical('------------------------------------------------')
+    log.critical('Statistics of output file:')
+    log.critical('  Name: %s' %swwname)
+    log.critical('  Reference:')
+    log.critical('    Lower left corner: [%f, %f]'
+                 % (geo_ref.get_xllcorner(), geo_ref.get_yllcorner()))
+    log.critical(' Start time: %f' %starttime)
+    log.critical('    Min time: %f' %mint)
+    log.critical('    Max time: %f' %maxt)
+    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 = outfile.variables[name][:]    # .flatten()
+        log.critical('    %s in [%f, %f]' % \
+                        (name, num.min(q), num.max(q)))        
+                 
+
+def _assert_lat_long(minlat, maxlat, minlon, maxlon):
+    """Check latitudes and longitudes for validity."""
+    
+    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
+
+    if minlat != None:
+        assert -90 < minlat < 90 , msg
+    if maxlat != None:
+        assert -90 < maxlat < 90 , msg
+        if minlat != None:
+            assert maxlat > minlat
+    if minlon != None:
+        assert -180 < minlon < 180 , msg
+    if maxlon != None:
+        assert -180 < maxlon < 180 , msg
+        if minlon != None:
+            assert maxlon > minlon
Index: /anuga_core/source/anuga/file_conversion/file_conversion.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/file_conversion.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/file_conversion.py	(revision 7742)
@@ -0,0 +1,966 @@
+""" Conversion routines.
+    ANUGA needs to deal with many different file formats, and this
+    module provides routines for easily converting between them.
+    
+    These routines are necessarily high level, sitting above the various
+    ANUGA modules. They take a file as input, and output a file.
+"""
+
+#non ANUGA imports
+from Scientific.IO.NetCDF import NetCDFFile
+import numpy as num
+import os.path
+import os
+
+#ANUGA imports
+from anuga.coordinate_transforms.geo_reference import Geo_reference, \
+     write_NetCDF_georeference, ensure_geo_reference
+from anuga.abstract_2d_finite_volumes.pmesh2domain import \
+     pmesh_to_domain_instance
+from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+
+
+#shallow water imports
+from anuga.shallow_water.sww_file import Read_sww, Write_sww
+from anuga.shallow_water.shallow_water_domain import Domain
+from anuga.shallow_water.shallow_water_domain import Domain
+
+
+def sww2obj(basefilename, size):
+    """ Convert netcdf based data output to obj
+
+        Convert SWW data to OBJ data.
+        basefilename Stem of filename, needs size and extension added.
+        size The number of lines to write.
+    """
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    # Get NetCDF
+    FN = create_filename('.', basefilename, 'sww', size)
+    log.critical('Reading from %s' % FN)
+    fid = NetCDFFile(FN, netcdf_mode_r)  #Open existing file for read
+
+    # Get the variables
+    x = fid.variables['x']
+    y = fid.variables['y']
+    z = fid.variables['elevation']
+    time = fid.variables['time']
+    stage = fid.variables['stage']
+
+    M = size  #Number of lines
+    xx = num.zeros((M,3), num.float)
+    yy = num.zeros((M,3), num.float)
+    zz = num.zeros((M,3), num.float)
+
+    for i in range(M):
+        for j in range(3):
+            xx[i,j] = x[i+j*M]
+            yy[i,j] = y[i+j*M]
+            zz[i,j] = z[i+j*M]
+
+    # Write obj for bathymetry
+    FN = create_filename('.', basefilename, 'obj', size)
+    write_obj(FN,xx,yy,zz)
+
+    # Now read all the data with variable information, combine with
+    # x,y info and store as obj
+    for k in range(len(time)):
+        t = time[k]
+        log.critical('Processing timestep %f' % t)
+
+        for i in range(M):
+            for j in range(3):
+                zz[i,j] = stage[k,i+j*M]
+
+        #Write obj for variable data
+        #FN = create_filename(basefilename, 'obj', size, time=t)
+        FN = create_filename('.', basefilename[:5], 'obj', size, time=t)
+        write_obj(FN, xx, yy, zz)
+
+
+##
+# @brief 
+# @param basefilename Stem of filename, needs size and extension added.
+def dat2obj(basefilename):
+    """Convert line based data output to obj
+    FIXME: Obsolete?
+    """
+
+    import glob, os
+    from anuga.config import data_dir
+
+    # Get bathymetry and x,y's
+    lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines()
+
+    M = len(lines)  #Number of lines
+    x = num.zeros((M,3), num.float)
+    y = num.zeros((M,3), num.float)
+    z = num.zeros((M,3), num.float)
+
+    for i, line in enumerate(lines):
+        tokens = line.split()
+        values = map(float, tokens)
+
+        for j in range(3):
+            x[i,j] = values[j*3]
+            y[i,j] = values[j*3+1]
+            z[i,j] = values[j*3+2]
+
+    # Write obj for bathymetry
+    write_obj(data_dir + os.sep + basefilename + '_geometry', x, y, z)
+
+    # Now read all the data files with variable information, combine with
+    # x,y info and store as obj.
+
+    files = glob.glob(data_dir + os.sep + basefilename + '*.dat')
+    for filename in files:
+        log.critical('Processing %s' % filename)
+
+        lines = open(data_dir + os.sep + filename, 'r').readlines()
+        assert len(lines) == M
+        root, ext = os.path.splitext(filename)
+
+        # Get time from filename
+        i0 = filename.find('_time=')
+        if i0 == -1:
+            #Skip bathymetry file
+            continue
+
+        i0 += 6  #Position where time starts
+        i1 = filename.find('.dat')
+
+        if i1 > i0:
+            t = float(filename[i0:i1])
+        else:
+            raise DataTimeError, 'Hmmmm'
+
+        for i, line in enumerate(lines):
+            tokens = line.split()
+            values = map(float,tokens)
+
+            for j in range(3):
+                z[i,j] = values[j]
+
+        # Write obj for variable data
+        write_obj(data_dir + os.sep + basefilename + '_time=%.4f' % t, x, y, z)
+
+##
+# @brief Convert time-series text file to TMS file.
+# @param filename 
+# @param quantity_names 
+# @param time_as_seconds 
+def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False):
+    """Template for converting typical text files with time series to
+    NetCDF tms file.
+
+    The file format is assumed to be either two fields separated by a comma:
+
+        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
+
+    E.g
+
+      31/08/04 00:00:00, 1.328223 0 0
+      31/08/04 00:15:00, 1.292912 0 0
+
+    or time (seconds), value0 value1 value2 ...
+
+      0.0, 1.328223 0 0
+      0.1, 1.292912 0 0
+
+    will provide a time dependent function f(t) with three attributes
+
+    filename is assumed to be the rootname with extenisons .txt and .sww
+    """
+
+    import time, calendar
+    from anuga.config import time_format
+    from anuga.utilities.numerical_tools import ensure_numeric
+
+    file_text = filename + '.txt'
+    fid = open(file_text)
+    line = fid.readline()
+    fid.close()
+
+    fields = line.split(',')
+    msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \
+          % file_text
+    assert len(fields) == 2, msg
+
+    if not time_as_seconds:
+        try:
+            starttime = calendar.timegm(time.strptime(fields[0], time_format))
+        except ValueError:
+            msg = 'First field in file %s must be' % file_text
+            msg += ' date-time with format %s.\n' % time_format
+            msg += 'I got %s instead.' % fields[0]
+            raise DataTimeError, msg
+    else:
+        try:
+            starttime = float(fields[0])
+        except Error:
+            msg = "Bad time format"
+            raise DataTimeError, msg
+
+    # Split values
+    values = []
+    for value in fields[1].split():
+        values.append(float(value))
+
+    q = ensure_numeric(values)
+
+    msg = 'ERROR: File must contain at least one independent value'
+    assert len(q.shape) == 1, msg
+
+    # Read times proper
+    from anuga.config import time_format
+    import time, calendar
+
+    fid = open(file_text)
+    lines = fid.readlines()
+    fid.close()
+
+    N = len(lines)
+    d = len(q)
+
+    T = num.zeros(N, num.float)       # Time
+    Q = num.zeros((N, d), num.float)  # Values
+
+    for i, line in enumerate(lines):
+        fields = line.split(',')
+        if not time_as_seconds:
+            realtime = calendar.timegm(time.strptime(fields[0], time_format))
+        else:
+            realtime = float(fields[0])
+        T[i] = realtime - starttime
+
+        for j, value in enumerate(fields[1].split()):
+            Q[i, j] = float(value)
+
+    msg = 'File %s must list time as a monotonuosly ' % filename
+    msg += 'increasing sequence'
+    assert num.alltrue(T[1:] - T[:-1] > 0), msg
+
+    #Create NetCDF file
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
+
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'Time series'
+
+    #Reference point
+    #Start time in seconds since the epoch (midnight 1/1/1970)
+    #FIXME: Use Georef
+    fid.starttime = starttime
+
+    # dimension definitions
+    #fid.createDimension('number_of_volumes', self.number_of_volumes)
+    #fid.createDimension('number_of_vertices', 3)
+
+    fid.createDimension('number_of_timesteps', len(T))
+
+    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
+
+    fid.variables['time'][:] = T
+
+    for i in range(Q.shape[1]):
+        try:
+            name = quantity_names[i]
+        except:
+            name = 'Attribute%d' % i
+
+        fid.createVariable(name, netcdf_float, ('number_of_timesteps',))
+        fid.variables[name][:] = Q[:,i]
+
+    fid.close()
+    
+    
+
+##
+# @brief 
+# @param filename
+# @param verbose
+def tsh2sww(filename, verbose=False):
+    """
+    to check if a tsh/msh file 'looks' good.
+    """
+
+    if verbose == True: log.critical('Creating domain from %s' % filename)
+
+    domain = pmesh_to_domain_instance(filename, Domain)
+
+    if verbose == True: log.critical("Number of triangles = %s" % len(domain))
+
+    domain.smooth = True
+    domain.format = 'sww'   #Native netcdf visualisation format
+    file_path, filename = path.split(filename)
+    filename, ext = path.splitext(filename)
+    domain.set_name(filename)
+    domain.reduction = mean
+
+    if verbose == True: log.critical("file_path = %s" % file_path)
+
+    if file_path == "":
+        file_path = "."
+    domain.set_datadir(file_path)
+
+    if verbose == True:
+        log.critical("Output written to %s%s%s.%s"
+                     % (domain.get_datadir(), sep, domain.get_name(),
+                        domain.format))
+
+    sww = SWW_file(domain)
+    sww.store_connectivity()
+    sww.store_timestep()
+
+
+##
+# @brief Convert CSIRO ESRI file to an SWW boundary file.
+# @param bath_dir 
+# @param elevation_dir 
+# @param ucur_dir 
+# @param vcur_dir 
+# @param sww_file 
+# @param minlat 
+# @param maxlat 
+# @param minlon 
+# @param maxlon 
+# @param zscale 
+# @param mean_stage 
+# @param fail_on_NaN 
+# @param elevation_NaN_filler 
+# @param bath_prefix 
+# @param elevation_prefix 
+# @param verbose 
+# @note Also convert latitude and longitude to UTM. All coordinates are
+#       assumed to be given in the GDA94 datum.
+def asc_csiro2sww(bath_dir,
+                  elevation_dir,
+                  ucur_dir,
+                  vcur_dir,
+                  sww_file,
+                  minlat=None, maxlat=None,
+                  minlon=None, maxlon=None,
+                  zscale=1,
+                  mean_stage=0,
+                  fail_on_NaN=True,
+                  elevation_NaN_filler=0,
+                  bath_prefix='ba',
+                  elevation_prefix='el',
+                  verbose=False):
+    """
+    Produce an sww boundary file, from esri ascii data from CSIRO.
+
+    Also convert latitude and longitude to UTM. All coordinates are
+    assumed to be given in the GDA94 datum.
+
+    assume:
+    All files are in esri ascii format
+
+    4 types of information
+    bathymetry
+    elevation
+    u velocity
+    v velocity
+
+    Assumptions
+    The metadata of all the files is the same
+    Each type is in a seperate directory
+    One bath file with extention .000
+    The time period is less than 24hrs and uniform.
+    """
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    from anuga.coordinate_transforms.redfearn import redfearn
+
+    # So if we want to change the precision it's done here
+    precision = netcdf_float 
+
+    # go in to the bath dir and load the only file,
+    bath_files = os.listdir(bath_dir)
+    bath_file = bath_files[0]
+    bath_dir_file =  bath_dir + os.sep + bath_file
+    bath_metadata, bath_grid =  _read_asc(bath_dir_file)
+
+    #Use the date.time of the bath file as a basis for
+    #the start time for other files
+    base_start = bath_file[-12:]
+
+    #go into the elevation dir and load the 000 file
+    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
+                         + base_start
+
+    elevation_files = os.listdir(elevation_dir)
+    ucur_files = os.listdir(ucur_dir)
+    vcur_files = os.listdir(vcur_dir)
+    elevation_files.sort()
+
+    # the first elevation file should be the
+    # file with the same base name as the bath data
+    assert elevation_files[0] == 'el' + base_start
+
+    number_of_latitudes = bath_grid.shape[0]
+    number_of_longitudes = bath_grid.shape[1]
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
+
+    longitudes = [bath_metadata['xllcorner'] + x*bath_metadata['cellsize'] \
+                  for x in range(number_of_longitudes)]
+    latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \
+                 for y in range(number_of_latitudes)]
+
+     # reverse order of lat, so the first lat represents the first grid row
+    latitudes.reverse()
+
+    kmin, kmax, lmin, lmax = get_min_max_indexes(latitudes[:],longitudes[:],
+                                                  minlat=minlat, maxlat=maxlat,
+                                                  minlon=minlon, maxlon=maxlon)
+
+    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
+    latitudes = latitudes[kmin:kmax]
+    longitudes = longitudes[lmin:lmax]
+    number_of_latitudes = len(latitudes)
+    number_of_longitudes = len(longitudes)
+    number_of_times = len(os.listdir(elevation_dir))
+    number_of_points = number_of_latitudes * number_of_longitudes
+    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
+
+    #Work out the times
+    if len(elevation_files) > 1:
+        # Assume: The time period is less than 24hrs.
+        time_period = (int(elevation_files[1][-3:]) \
+                       - int(elevation_files[0][-3:])) * 60*60
+        times = [x*time_period for x in range(len(elevation_files))]
+    else:
+        times = [0.0]
+
+    if verbose:
+        log.critical('------------------------------------------------')
+        log.critical('Statistics:')
+        log.critical('  Extent (lat/lon):')
+        log.critical('    lat in [%f, %f], len(lat) == %d'
+                     % (min(latitudes), max(latitudes), len(latitudes)))
+        log.critical('    lon in [%f, %f], len(lon) == %d'
+                     % (min(longitudes), max(longitudes), len(longitudes)))
+        log.critical('    t in [%f, %f], len(t) == %d'
+                     % (min(times), max(times), len(times)))
+
+    ######### WRITE THE SWW FILE #############
+
+    # NetCDF file definition
+    outfile = NetCDFFile(sww_file, netcdf_mode_w)
+
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'Converted from XXX'
+
+    #For sww compatibility
+    outfile.smoothing = 'Yes'
+    outfile.order = 1
+
+    #Start time in seconds since the epoch (midnight 1/1/1970)
+    outfile.starttime = starttime = times[0]
+
+    # dimension definitions
+    outfile.createDimension('number_of_volumes', number_of_volumes)
+    outfile.createDimension('number_of_vertices', 3)
+    outfile.createDimension('number_of_points', number_of_points)
+    outfile.createDimension('number_of_timesteps', number_of_times)
+
+    # variable definitions
+    outfile.createVariable('x', precision, ('number_of_points',))
+    outfile.createVariable('y', precision, ('number_of_points',))
+    outfile.createVariable('elevation', precision, ('number_of_points',))
+
+    #FIXME: Backwards compatibility
+    #outfile.createVariable('z', precision, ('number_of_points',))
+    #################################
+
+    outfile.createVariable('volumes', netcdf_int, ('number_of_volumes',
+                                                   'number_of_vertices'))
+
+    outfile.createVariable('time', precision, ('number_of_timesteps',))
+
+    outfile.createVariable('stage', precision, ('number_of_timesteps',
+                                                'number_of_points'))
+
+    outfile.createVariable('xmomentum', precision, ('number_of_timesteps',
+                                                    'number_of_points'))
+
+    outfile.createVariable('ymomentum', precision, ('number_of_timesteps',
+                                                    'number_of_points'))
+
+    #Store
+    from anuga.coordinate_transforms.redfearn import redfearn
+
+    x = num.zeros(number_of_points, num.float)  #Easting
+    y = num.zeros(number_of_points, num.float)  #Northing
+
+    if verbose: log.critical('Making triangular grid')
+
+    #Get zone of 1st point.
+    refzone, _, _ = redfearn(latitudes[0], longitudes[0])
+
+    vertices = {}
+    i = 0
+    for k, lat in enumerate(latitudes):
+        for l, lon in enumerate(longitudes):
+            vertices[l,k] = i
+
+            zone, easting, northing = redfearn(lat, lon)
+
+            #msg = 'Zone boundary crossed at longitude =', lon
+            #assert zone == refzone, msg
+            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
+            x[i] = easting
+            y[i] = northing
+            i += 1
+
+    #Construct 2 triangles per 'rectangular' element
+    volumes = []
+    for l in range(number_of_longitudes-1):    #X direction
+        for k in range(number_of_latitudes-1): #Y direction
+            v1 = vertices[l,k+1]
+            v2 = vertices[l,k]
+            v3 = vertices[l+1,k+1]
+            v4 = vertices[l+1,k]
+
+            #Note, this is different to the ferrit2sww code
+            #since the order of the lats is reversed.
+            volumes.append([v1,v3,v2]) #Upper element
+            volumes.append([v4,v2,v3]) #Lower element
+
+    volumes = num.array(volumes, num.int)      #array default#
+
+    geo_ref = Geo_reference(refzone, min(x), min(y))
+    geo_ref.write_NetCDF(outfile)
+
+    # This will put the geo ref in the middle
+    #geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.)
+
+    if verbose:
+        log.critical('------------------------------------------------')
+        log.critical('More Statistics:')
+        log.critical('  Extent (/lon):')
+        log.critical('    x in [%f, %f], len(lat) == %d'
+                     % (min(x), max(x), len(x)))
+        log.critical('    y in [%f, %f], len(lon) == %d'
+                     % (min(y), max(y), len(y)))
+        log.critical('geo_ref: ', geo_ref)
+
+    z = num.resize(bath_grid,outfile.variables['elevation'][:].shape)
+    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
+    # FIXME (Ole): Remove once viewer has been recompiled and changed
+    #              to use elevation instead of z
+    #outfile.variables['z'][:] = z
+    outfile.variables['elevation'][:] = z
+    outfile.variables['volumes'][:] = volumes.astype(num.int32) # On Opteron 64
+
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum = outfile.variables['ymomentum']
+
+    outfile.variables['time'][:] = times   #Store time relative
+
+    if verbose: log.critical('Converting quantities')
+
+    n = number_of_times
+    for j in range(number_of_times):
+        # load in files
+        elevation_meta, elevation_grid = \
+            _read_asc(elevation_dir + os.sep + elevation_files[j])
+
+        _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j])
+        _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j])
+
+        #cut matrix to desired size
+        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
+        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
+        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
+
+        # handle missing values
+        missing = (elevation_grid == elevation_meta['NODATA_value'])
+        if num.sometrue (missing):
+            if fail_on_NaN:
+                msg = 'File %s contains missing values' \
+                      % (elevation_files[j])
+                raise DataMissingValuesError, msg
+            else:
+                elevation_grid = elevation_grid*(missing==0) \
+                                 + missing*elevation_NaN_filler
+
+        if verbose and j % ((n+10)/10) == 0: log.critical('  Doing %d of %d'
+                                                          % (j, n))
+
+        i = 0
+        for k in range(number_of_latitudes):      #Y direction
+            for l in range(number_of_longitudes): #X direction
+                w = zscale*elevation_grid[k,l] + mean_stage
+                stage[j,i] = w
+                h = w - z[i]
+                xmomentum[j,i] = u_momentum_grid[k,l]*h
+                ymomentum[j,i] = v_momentum_grid[k,l]*h
+                i += 1
+
+    outfile.close()
+
+
+
+def _read_asc(filename, verbose=False):
+    """Read esri file from the following ASCII format (.asc)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    filename Path to the file to read.
+    verbose True if this function is to be verbose.
+    """
+
+    datafile = open(filename)
+
+    if verbose: log.critical('Reading DEM from %s' % filename)
+
+    lines = datafile.readlines()
+    datafile.close()
+
+    if verbose: log.critical('Got %d lines' % len(lines))
+
+    ncols = int(lines.pop(0).split()[1].strip())
+    nrows = int(lines.pop(0).split()[1].strip())
+    xllcorner = float(lines.pop(0).split()[1].strip())
+    yllcorner = float(lines.pop(0).split()[1].strip())
+    cellsize = float(lines.pop(0).split()[1].strip())
+    NODATA_value = float(lines.pop(0).split()[1].strip())
+
+    assert len(lines) == nrows
+
+    #Store data
+    grid = []
+
+    n = len(lines)
+    for i, line in enumerate(lines):
+        cells = line.split()
+        assert len(cells) == ncols
+        grid.append(num.array([float(x) for x in cells]))
+    grid = num.array(grid)
+
+    return {'xllcorner':xllcorner,
+            'yllcorner':yllcorner,
+            'cellsize':cellsize,
+            'NODATA_value':NODATA_value}, grid
+
+
+##
+# @brief Convert URS file to SWW file.
+# @param basename_in Stem of the input filename.
+# @param basename_out Stem of the output filename.
+# @param verbose True if this function is to be verbose.
+# @param remove_nc_files 
+# @param minlat Sets extent of area to be used.  If not supplied, full extent.
+# @param maxlat Sets extent of area to be used.  If not supplied, full extent.
+# @param minlon Sets extent of area to be used.  If not supplied, full extent.
+# @param maxlon Sets extent of area to be used.  If not supplied, full extent.
+# @param mint 
+# @param maxt 
+# @param mean_stage 
+# @param origin A 3-tuple with geo referenced UTM coordinates
+# @param zscale 
+# @param fail_on_NaN 
+# @param NaN_filler 
+# @param elevation 
+# @note Also convert latitude and longitude to UTM. All coordinates are
+#       assumed to be given in the GDA94 datum.
+def urs2sww(basename_in='o', basename_out=None, verbose=False,
+            remove_nc_files=True,
+            minlat=None, maxlat=None,
+            minlon=None, maxlon=None,
+            mint=None, maxt=None,
+            mean_stage=0,
+            origin=None,
+            zscale=1,
+            fail_on_NaN=True,
+            NaN_filler=0):
+    """Convert a URS file to an SWW file.
+    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-mux2, basefilename-e-mux2 and
+    basefilename-n-mux2 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 for  a grid.
+
+    min's and max's: If omitted - full extend is used.
+    To include a value min may equal it, while max must exceed it.
+    Lat and lon are assumed to be in decimal degrees.
+    NOTE: minlon is the most east boundary.
+
+    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.
+
+    URS C binary format has data orgainised 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.
+    """
+
+    if basename_out == None:
+        basename_out = basename_in
+
+    files_out = urs2nc(basename_in, basename_out)
+
+    ferret2sww(basename_out,
+               minlat=minlat,
+               maxlat=maxlat,
+               minlon=minlon,
+               maxlon=maxlon,
+               mint=mint,
+               maxt=maxt,
+               mean_stage=mean_stage,
+               origin=origin,
+               zscale=zscale,
+               fail_on_NaN=fail_on_NaN,
+               NaN_filler=NaN_filler,
+               inverted_bathymetry=True,
+               verbose=verbose)
+    
+    if remove_nc_files:
+        for file_out in files_out:
+            os.remove(file_out)
+
+
+##
+# @brief Convert 3 URS files back to 4 NC files.
+# @param basename_in Stem of the input filenames.
+# @param basename_outStem of the output filenames.
+# @note The name of the urs file names must be:
+#          [basename_in]-z-mux
+#          [basename_in]-e-mux
+#          [basename_in]-n-mux
+def urs2nc(basename_in='o', basename_out='urs'):
+    """Convert the 3 urs files to 4 nc files.
+
+    The name of the urs file names must be;
+    [basename_in]-z-mux
+    [basename_in]-e-mux
+    [basename_in]-n-mux
+    """
+
+    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
+                basename_in + EAST_VELOCITY_LABEL,
+                basename_in + NORTH_VELOCITY_LABEL]
+    files_out = [basename_out + '_ha.nc',
+                 basename_out + '_ua.nc',
+                 basename_out + '_va.nc']
+    quantities = ['HA', 'UA', 'VA']
+
+    #if os.access(files_in[0]+'.mux', os.F_OK) == 0 :
+    for i, file_name in enumerate(files_in):
+        if os.access(file_name, os.F_OK) == 0:
+            if os.access(file_name + '.mux', os.F_OK) == 0 :
+                msg = 'File %s does not exist or is not accessible' % file_name
+                raise IOError, msg
+            else:
+               files_in[i] += '.mux'
+               log.critical("file_name %s" % file_name)
+
+    hashed_elevation = None
+    for file_in, file_out, quantity in map(None, files_in,
+                                           files_out,
+                                           quantities):
+        lonlatdep, lon, lat, depth = _binary_c2nc(file_in,
+                                                  file_out,
+                                                  quantity)
+        if hashed_elevation == None:
+            elevation_file = basename_out + '_e.nc'
+            write_elevation_nc(elevation_file,
+                               lon,
+                               lat,
+                               depth)
+            hashed_elevation = myhash(lonlatdep)
+        else:
+            msg = "The elevation information in the mux files is inconsistent"
+            assert hashed_elevation == myhash(lonlatdep), msg
+
+    files_out.append(elevation_file)
+
+    return files_out
+
+
+##
+# @brief Convert a quantity URS file to a NetCDF file.
+# @param file_in Path to input URS file.
+# @param file_out Path to the output file.
+# @param quantity Name of the quantity to be written to the output file.
+# @return A tuple (lonlatdep, lon, lat, depth).
+def _binary_c2nc(file_in, file_out, quantity):
+    """Reads in a quantity urs file and writes a quantity nc file.
+    Additionally, returns the depth and lat, long info,
+    so it can be written to a file.
+    """
+
+    columns = 3                           # long, lat , depth
+    mux_file = open(file_in, 'rb')
+
+    # Number of points/stations
+    (points_num,) = unpack('i', mux_file.read(4))
+
+    # nt, int - Number of time steps
+    (time_step_count,) = unpack('i', mux_file.read(4))
+
+    #dt, float - time step, seconds
+    (time_step,) = unpack('f', mux_file.read(4))
+
+    msg = "Bad data in the mux file."
+    if points_num < 0:
+        mux_file.close()
+        raise ANUGAError, msg
+    if time_step_count < 0:
+        mux_file.close()
+        raise ANUGAError, msg
+    if time_step < 0:
+        mux_file.close()
+        raise ANUGAError, msg
+
+    lonlatdep = p_array.array('f')
+    lonlatdep.read(mux_file, columns * points_num)
+    lonlatdep = num.array(lonlatdep, dtype=num.float)
+    lonlatdep = num.reshape(lonlatdep, (points_num, columns))
+
+    lon, lat, depth = lon_lat2grid(lonlatdep)
+    lon_sorted = list(lon)
+    lon_sorted.sort()
+
+    if not num.alltrue(lon == lon_sorted):
+        msg = "Longitudes in mux file are not in ascending order"
+        raise IOError, msg
+
+    lat_sorted = list(lat)
+    lat_sorted.sort()
+
+    nc_file = Write_nc(quantity,
+                       file_out,
+                       time_step_count,
+                       time_step,
+                       lon,
+                       lat)
+
+    for i in range(time_step_count):
+        #Read in a time slice from mux file
+        hz_p_array = p_array.array('f')
+        hz_p_array.read(mux_file, points_num)
+        hz_p = num.array(hz_p_array, dtype=num.float)
+        hz_p = num.reshape(hz_p, (len(lon), len(lat)))
+        hz_p = num.transpose(hz_p)  # mux has lat varying fastest, nc has long v.f.
+
+        #write time slice to nc file
+        nc_file.store_timestep(hz_p)
+
+    mux_file.close()
+    nc_file.close()
+
+    return lonlatdep, lon, lat, depth
+
+
+
+##
+# @brief Return max&min indexes (for slicing) of area specified.
+# @param latitudes_ref ??
+# @param longitudes_ref ??
+# @param minlat Minimum latitude of specified area.
+# @param maxlat Maximum latitude of specified area.
+# @param minlon Minimum longitude of specified area.
+# @param maxlon Maximum longitude of specified area.
+# @return Tuple (lat_min_index, lat_max_index, lon_min_index, lon_max_index)
+def get_min_max_indexes(latitudes_ref, longitudes_ref,
+                         minlat=None, maxlat=None,
+                         minlon=None, maxlon=None):
+    """
+    Return max, min indexes (for slicing) of the lat and long arrays to cover
+    the area specified with min/max lat/long.
+
+    Think of the latitudes and longitudes describing a 2d surface.
+    The area returned is, if possible, just big enough to cover the
+    inputed max/min area. (This will not be possible if the max/min area
+    has a section outside of the latitudes/longitudes area.)
+
+    asset  longitudes are sorted,
+    long - from low to high (west to east, eg 148 - 151)
+    assert latitudes are sorted, ascending or decending
+    """
+
+    latitudes = latitudes_ref[:]
+    longitudes = longitudes_ref[:]
+
+    latitudes = ensure_numeric(latitudes)
+    longitudes = ensure_numeric(longitudes)
+
+    assert num.allclose(num.sort(longitudes), longitudes)
+
+    #print latitudes[0],longitudes[0]
+    #print len(latitudes),len(longitudes)
+    #print latitudes[len(latitudes)-1],longitudes[len(longitudes)-1]
+
+    lat_ascending = True
+    if not num.allclose(num.sort(latitudes), latitudes):
+        lat_ascending = False
+        # reverse order of lat, so it's in ascending order
+        latitudes = latitudes[::-1]
+        assert num.allclose(num.sort(latitudes), latitudes)
+
+    largest_lat_index = len(latitudes)-1
+
+    #Cut out a smaller extent.
+    if minlat == None:
+        lat_min_index = 0
+    else:
+        lat_min_index = num.searchsorted(latitudes, minlat)-1
+        if lat_min_index < 0:
+            lat_min_index = 0
+
+    if maxlat == None:
+        lat_max_index = largest_lat_index #len(latitudes)
+    else:
+        lat_max_index = num.searchsorted(latitudes, maxlat)
+        if lat_max_index > largest_lat_index:
+            lat_max_index = largest_lat_index
+
+    if minlon == None:
+        lon_min_index = 0
+    else:
+        lon_min_index = num.searchsorted(longitudes, minlon)-1
+        if lon_min_index < 0:
+            lon_min_index = 0
+
+    if maxlon == None:
+        lon_max_index = len(longitudes)
+    else:
+        lon_max_index = num.searchsorted(longitudes, maxlon)
+
+    # Reversing the indexes, if the lat array is decending
+    if lat_ascending is False:
+        lat_min_index, lat_max_index = largest_lat_index - lat_max_index, \
+                                       largest_lat_index - lat_min_index
+    lat_max_index = lat_max_index + 1 # taking into account how slicing works
+    lon_max_index = lon_max_index + 1 # taking into account how slicing works
+
+    return lat_min_index, lat_max_index, lon_min_index, lon_max_index
+
+
Index: /anuga_core/source/anuga/file_conversion/test_all.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/test_all.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/test_all.py	(revision 7742)
@@ -0,0 +1,83 @@
+"""Regression testing framework
+This module will search for scripts in the same directory named
+test_*.py.  Each such script should be a test suite that tests a
+module through PyUnit. This script will aggregate all
+found test suites into one big test suite and run them all at once.
+"""
+
+# Author: Mark Pilgrim
+# Modified by Ole Nielsen
+
+import unittest
+import os
+
+
+#List files that should be excluded from the testing process.
+#E.g. if they are known to fail and under development
+exclude = [] 
+
+
+def get_test_files(path):
+
+    import sys
+
+    files = os.listdir(path)
+
+    #Check sub directories
+    test_files = []
+    for file in files:
+        if os.path.isdir(file):
+            sys.path.append(file)
+            #print 'Recursing into', file
+            test_files += get_test_files(path + os.sep + file)
+        elif file[:5] == 'test_' and file[-2:] == 'py':
+            #print 'Appending', file
+            test_files.append(file)
+        else:
+            pass
+    return test_files
+
+
+
+def regressionTest():
+    import sys, os, re, unittest
+    path = os.path.split(sys.argv[0])[0] or os.getcwd()
+
+
+    files = get_test_files(path)
+
+
+
+    #test = re.compile('^test_[\w]*.py$', re.IGNORECASE)
+    #files = filter(test.search, files)
+
+
+    try:
+        files.remove(__file__)  #Remove self from list (Ver 2.3. or later)
+    except:
+        files.remove('test_all.py')
+
+    print 'Testing:'
+    for file in files:
+        print '  ' + file
+
+    if globals().has_key('exclude'):
+        for file in exclude:
+            files.remove(file)
+            print 'WARNING: File '+ file + ' excluded from testing'
+
+
+    filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    #print "files",files 
+    moduleNames = map(filenameToModuleName, files)
+    #print "moduleNames",moduleNames
+    modules = map(__import__, moduleNames)
+    load = unittest.defaultTestLoader.loadTestsFromModule
+    return unittest.TestSuite(map(load, modules))
+
+################################################################################
+
+if __name__ == '__main__':   
+    suite = regressionTest()
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
Index: /anuga_core/source/anuga/file_conversion/test_dem2pts.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/test_dem2pts.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/test_dem2pts.py	(revision 7742)
@@ -0,0 +1,399 @@
+import sys
+import unittest
+import numpy as num
+import copy
+import os
+
+# ANUGA modules
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
+                            netcdf_float
+
+from dem2pts import dem2pts
+from asc2dem import asc2dem
+
+class Test_Dem2Pts(unittest.TestCase):
+    """ A suite of tests to test file conversion functions.
+        These tests are quite coarse-grained: converting a file
+        and checking that its headers and some of its contents
+        are correct.
+    """ 
+ 
+    def test_dem2pts_bounding_box_v2(self):
+        """Test conversion from dem in ascii format to native NetCDF format
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test asc file
+        root = 'demtest'
+
+        filename = root+'.asc'
+        fid = open(filename, 'w')
+        fid.write("""ncols         10
+nrows         10
+xllcorner     2000
+yllcorner     3000
+cellsize      1
+NODATA_value  -9999
+""")
+        #Create linear function
+        ref_points = []
+        ref_elevation = []
+        x0 = 2000
+        y = 3010
+        yvec = range(10)
+        xvec = range(10)
+        z = -1
+        for i in range(10):
+            y = y - 1
+            for j in range(10):
+                x = x0 + xvec[j]
+                z += 1
+                ref_points.append ([x,y])
+                ref_elevation.append(z)
+                fid.write('%f ' %z)
+            fid.write('\n')
+
+        fid.close()
+
+        #print 'sending pts', ref_points
+        #print 'sending elev', ref_elevation
+
+        #Write prj file with metadata
+        metafilename = root+'.prj'
+        fid = open(metafilename, 'w')
+
+
+        fid.write("""Projection UTM
+Zone 56
+Datum WGS84
+Zunits NO
+Units METERS
+Spheroid WGS84
+Xshift 0.0000000000
+Yshift 10000000.0000000000
+Parameters
+""")
+        fid.close()
+
+        #Convert to NetCDF pts
+        asc2dem(root)
+        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
+                northing_min=3003.0, northing_max=3006.0,
+                verbose=False)
+
+        #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
+        assert fid.xllcorner[0] == 2002.0
+        assert fid.yllcorner[0] == 3003.0
+
+        #create new reference points
+        newz = []
+        newz[0:5] = ref_elevation[32:38]
+        newz[6:11] = ref_elevation[42:48]
+        newz[12:17] = ref_elevation[52:58]
+        newz[18:23] = ref_elevation[62:68]
+        ref_elevation = []
+        ref_elevation = newz
+        ref_points = []
+        x0 = 2002
+        y = 3007
+        yvec = range(4)
+        xvec = range(6)
+        for i in range(4):
+            y = y - 1
+            ynew = y - 3003.0
+            for j in range(6):
+                x = x0 + xvec[j]
+                xnew = x - 2002.0
+                ref_points.append ([xnew,ynew]) #Relative point values
+
+        assert num.allclose(points, ref_points)
+
+        assert num.allclose(elevation, ref_elevation)
+
+        #Cleanup
+        fid.close()
+
+
+        os.remove(root + '.pts')
+        os.remove(root + '.dem')
+        os.remove(root + '.asc')
+        os.remove(root + '.prj')
+
+
+    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
+        """Test conversion from dem in ascii format to native NetCDF format
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test asc file
+        root = 'demtest'
+
+        filename = root+'.asc'
+        fid = open(filename, 'w')
+        fid.write("""ncols         10
+nrows         10
+xllcorner     2000
+yllcorner     3000
+cellsize      1
+NODATA_value  -9999
+""")
+        #Create linear function
+        ref_points = []
+        ref_elevation = []
+        x0 = 2000
+        y = 3010
+        yvec = range(10)
+        xvec = range(10)
+        #z = range(100)
+        z = num.zeros(100, num.int)     #array default#
+        NODATA_value = -9999
+        count = -1
+        for i in range(10):
+            y = y - 1
+            for j in range(10):
+                x = x0 + xvec[j]
+                ref_points.append ([x,y])
+                count += 1
+                z[count] = (4*i - 3*j)%13
+                if j == 4: z[count] = NODATA_value #column inside clipping region
+                if j == 8: z[count] = NODATA_value #column outside clipping region
+                if i == 9: z[count] = NODATA_value #row outside clipping region
+                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
+                ref_elevation.append( z[count] )
+                fid.write('%f ' %z[count])
+            fid.write('\n')
+
+        fid.close()
+
+        #print 'sending elev', ref_elevation
+
+        #Write prj file with metadata
+        metafilename = root+'.prj'
+        fid = open(metafilename, 'w')
+
+
+        fid.write("""Projection UTM
+Zone 56
+Datum WGS84
+Zunits NO
+Units METERS
+Spheroid WGS84
+Xshift 0.0000000000
+Yshift 10000000.0000000000
+Parameters
+""")
+        fid.close()
+
+        #Convert to NetCDF pts
+        asc2dem(root)
+        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
+                northing_min=3003.0, northing_max=3006.0)
+
+        #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
+        assert fid.xllcorner[0] == 2002.0
+        assert fid.yllcorner[0] == 3003.0
+
+        #create new reference points
+        newz = num.zeros(19, num.int)       #array default#
+        newz[0:2] = ref_elevation[32:34]
+        newz[2:5] = ref_elevation[35:38]
+        newz[5:7] = ref_elevation[42:44]
+        newz[7] = ref_elevation[45]
+        newz[8] = ref_elevation[47]
+        newz[9:11] = ref_elevation[52:54]
+        newz[11:14] = ref_elevation[55:58]
+        newz[14:16] = ref_elevation[62:64]
+        newz[16:19] = ref_elevation[65:68]
+
+
+        ref_elevation = newz
+        ref_points = []
+        new_ref_points = []
+        x0 = 2002
+        y = 3007
+        yvec = range(4)
+        xvec = range(6)
+        for i in range(4):
+            y = y - 1
+            ynew = y - 3003.0
+            for j in range(6):
+                x = x0 + xvec[j]
+                xnew = x - 2002.0
+                if j <> 2 and (i<>1 or j<>4):
+                    ref_points.append([x,y])
+                    new_ref_points.append ([xnew,ynew])
+
+
+        assert num.allclose(points, new_ref_points)
+        assert num.allclose(elevation, ref_elevation)
+
+        #Cleanup
+        fid.close()
+
+
+        os.remove(root + '.pts')
+        os.remove(root + '.dem')
+        os.remove(root + '.asc')
+        os.remove(root + '.prj')
+
+
+    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
+        """Test conversion from dem in ascii format to native NetCDF format
+        Check missing values on clipping boundary
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test asc file
+        root = 'demtest'
+
+        filename = root+'.asc'
+        fid = open(filename, 'w')
+        fid.write("""ncols         10
+nrows         10
+xllcorner     2000
+yllcorner     3000
+cellsize      1
+NODATA_value  -9999
+""")
+        #Create linear function
+        ref_points = []
+        ref_elevation = []
+        x0 = 2000
+        y = 3010
+        yvec = range(10)
+        xvec = range(10)
+        #z = range(100)
+        z = num.zeros(100, num.int)     #array default#
+        NODATA_value = -9999
+        count = -1
+        for i in range(10):
+            y = y - 1
+            for j in range(10):
+                x = x0 + xvec[j]
+                ref_points.append ([x,y])
+                count += 1
+                z[count] = (4*i - 3*j)%13
+                if j == 4: z[count] = NODATA_value #column inside clipping region
+                if j == 8: z[count] = NODATA_value #column outside clipping region
+                if i == 6: z[count] = NODATA_value #row on clipping boundary
+                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
+                ref_elevation.append( z[count] )
+                fid.write('%f ' %z[count])
+            fid.write('\n')
+
+        fid.close()
+
+        #print 'sending elev', ref_elevation
+
+        #Write prj file with metadata
+        metafilename = root+'.prj'
+        fid = open(metafilename, 'w')
+
+
+        fid.write("""Projection UTM
+Zone 56
+Datum WGS84
+Zunits NO
+Units METERS
+Spheroid WGS84
+Xshift 0.0000000000
+Yshift 10000000.0000000000
+Parameters
+""")
+        fid.close()
+
+        #Convert to NetCDF pts
+        asc2dem(root)
+        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
+                northing_min=3003.0, northing_max=3006.0)
+
+        #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
+        assert fid.xllcorner[0] == 2002.0
+        assert fid.yllcorner[0] == 3003.0
+
+        #create new reference points
+        newz = num.zeros(14, num.int)       #array default#
+        newz[0:2] = ref_elevation[32:34]
+        newz[2:5] = ref_elevation[35:38]
+        newz[5:7] = ref_elevation[42:44]
+        newz[7] = ref_elevation[45]
+        newz[8] = ref_elevation[47]
+        newz[9:11] = ref_elevation[52:54]
+        newz[11:14] = ref_elevation[55:58]
+
+
+
+        ref_elevation = newz
+        ref_points = []
+        new_ref_points = []
+        x0 = 2002
+        y = 3007
+        yvec = range(4)
+        xvec = range(6)
+        for i in range(4):
+            y = y - 1
+            ynew = y - 3003.0
+            for j in range(6):
+                x = x0 + xvec[j]
+                xnew = x - 2002.0
+                if j <> 2 and (i<>1 or j<>4) and i<>3:
+                    ref_points.append([x,y])
+                    new_ref_points.append ([xnew,ynew])
+
+
+        #print points[:],points[:].shape
+        #print new_ref_points, len(new_ref_points)
+
+        assert num.allclose(elevation, ref_elevation)
+        assert num.allclose(points, new_ref_points)
+
+
+        #Cleanup
+        fid.close()
+
+
+        os.remove(root + '.pts')
+        os.remove(root + '.dem')
+        os.remove(root + '.asc')
+        os.remove(root + '.prj')
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Dem2Pts,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
Index: /anuga_core/source/anuga/file_conversion/test_file_conversion.py
===================================================================
--- /anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7742)
+++ /anuga_core/source/anuga/file_conversion/test_file_conversion.py	(revision 7742)
@@ -0,0 +1,2249 @@
+
+from anuga.shallow_water.shallow_water_domain import Domain
+from ferret2sww import ferret2sww
+from anuga.utilities.numerical_tools import ensure_numeric, mean
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary
+from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+from anuga.shallow_water.sww_file import SWW_file
+from data_manager import extent_sww
+from anuga.config import netcdf_float, epsilon, g
+from Scientific.IO.NetCDF import NetCDFFile
+import sys
+import unittest
+import numpy as num
+import copy
+import os
+
+
+class Test_File_Conversion(unittest.TestCase):
+    """ A suite of tests to test file conversion functions.
+        These tests are quite coarse-grained: converting a file
+        and checking that its headers and some of its contents
+        are correct.
+    """
+    verbose = False
+
+    def set_verbose(self):
+        Test_File_Conversion.verbose = True
+        
+    def setUp(self):
+        import time
+        
+        self.verbose = Test_File_Conversion.verbose
+        # Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        # Set some field values
+        domain.set_quantity('elevation', lambda x,y: -x)
+        domain.set_quantity('friction', 0.03)
+
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        ######################
+        #Initial condition - with jumps
+        bed = domain.quantities['elevation'].vertex_values
+        stage = num.zeros(bed.shape, num.float)
+
+        h = 0.3
+        for i in range(stage.shape[0]):
+            if i % 2 == 0:
+                stage[i,:] = bed[i,:] + h
+            else:
+                stage[i,:] = bed[i,:]
+
+        domain.set_quantity('stage', stage)
+
+
+        domain.distribute_to_vertices_and_edges()               
+        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
+
+
+        self.domain = domain
+
+        C = domain.get_vertex_coordinates()
+        self.X = C[:,0:6:2].copy()
+        self.Y = C[:,1:6:2].copy()
+
+        self.F = bed
+
+        #Write A testfile (not realistic. Values aren't realistic)
+        self.test_MOST_file = 'most_small'
+
+        longitudes = [150.66667, 150.83334, 151., 151.16667]
+        latitudes = [-34.5, -34.33333, -34.16667, -34]
+
+        long_name = 'LON'
+        lat_name = 'LAT'
+
+        nx = 4
+        ny = 4
+        six = 6
+
+
+        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
+            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)
+
+            fid.createDimension(long_name,nx)
+            fid.createVariable(long_name,netcdf_float,(long_name,))
+            fid.variables[long_name].point_spacing='uneven'
+            fid.variables[long_name].units='degrees_east'
+            fid.variables[long_name].assignValue(longitudes)
+
+            fid.createDimension(lat_name,ny)
+            fid.createVariable(lat_name,netcdf_float,(lat_name,))
+            fid.variables[lat_name].point_spacing='uneven'
+            fid.variables[lat_name].units='degrees_north'
+            fid.variables[lat_name].assignValue(latitudes)
+
+            fid.createDimension('TIME',six)
+            fid.createVariable('TIME',netcdf_float,('TIME',))
+            fid.variables['TIME'].point_spacing='uneven'
+            fid.variables['TIME'].units='seconds'
+            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
+
+
+            name = ext[1:3].upper()
+            if name == 'E.': name = 'ELEVATION'
+            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
+            fid.variables[name].units='CENTIMETERS'
+            fid.variables[name].missing_value=-1.e+034
+
+            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
+                                              [-0.1214216, 0, 0, 0],
+                                              [0, 0, 0, 0],
+                                              [0, 0, 0, 0]],
+                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
+                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
+                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
+                                              [0, 2.291054e-005, 2.291054e-005, 0]],
+                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
+                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
+                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
+                                              [0, 0.0001374632, 0.0001374632, 0]],
+                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
+                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
+                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
+                                              [0, 0.0002520159, 0.0002520159, 0]],
+                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
+                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
+                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
+                                              [0, 0.0003665686, 0.0003665686, 0]],
+                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
+                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
+                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
+                                              [0, 0.0004811212, 0.0004811212, 0]]])
+
+
+            fid.close()
+
+
+
+
+    def tearDown(self):
+        import os
+        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
+            #print 'Trying to remove', self.test_MOST_file + ext
+            os.remove(self.test_MOST_file + ext)
+
+    def test_ferret2sww1(self):
+        """Test that georeferencing etc works when converting from
+        ferret format (lat/lon) to sww format (UTM)
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+        import os, sys
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
+        #
+        # First value (index=0) in small_ha.nc is 0.3400644 cm,
+        # Fourth value (index==3) is -6.50198 cm
+
+
+
+        #Read
+        from anuga.coordinate_transforms.redfearn import redfearn
+        #fid = NetCDFFile(self.test_MOST_file)
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+        first_value = fid.variables['HA'][:][0,0,0]
+        fourth_value = fid.variables['HA'][:][0,0,3]
+        fid.close()
+
+
+        #Call conversion (with zero origin)
+        #ferret2sww('small', verbose=False,
+        #           origin = (56, 0, 0))
+        ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 0, 0))
+
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)
+        #print zone, e, n
+
+        #Read output file 'small.sww'
+        #fid = NetCDFFile('small.sww')
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        #Check that first coordinate is correctly represented
+        assert num.allclose(x[0], e)
+        assert num.allclose(y[0], n)
+
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
+
+        assert num.allclose(stage[0,0], first_value/100)  #Meters
+
+        #Check fourth value
+        assert num.allclose(stage[0,3], fourth_value/100)  #Meters
+
+        fid.close()
+
+        #Cleanup
+        import os
+        os.remove(self.test_MOST_file + '.sww')
+
+
+
+    def test_ferret2sww_zscale(self):
+        """Test that zscale workse
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+        import os, sys
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
+        #
+        # First value (index=0) in small_ha.nc is 0.3400644 cm,
+        # Fourth value (index==3) is -6.50198 cm
+
+
+        #Read
+        from anuga.coordinate_transforms.redfearn import redfearn
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+        first_value = fid.variables['HA'][:][0,0,0]
+        fourth_value = fid.variables['HA'][:][0,0,3]
+        fid.close()
+
+        #Call conversion (with no scaling)
+        ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 0, 0))
+
+        #Work out the UTM coordinates for first point
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        #Check values
+        stage_1 = fid.variables['stage'][:]
+        xmomentum_1 = fid.variables['xmomentum'][:]
+        ymomentum_1 = fid.variables['ymomentum'][:]
+
+        assert num.allclose(stage_1[0,0], first_value/100)  #Meters
+        assert num.allclose(stage_1[0,3], fourth_value/100)  #Meters
+
+        fid.close()
+
+        #Call conversion (with scaling)
+        ferret2sww(self.test_MOST_file,
+                   zscale = 5,
+                   verbose=self.verbose,
+                   origin = (56, 0, 0))
+
+        #Work out the UTM coordinates for first point
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        #Check values
+        stage_5 = fid.variables['stage'][:]
+        xmomentum_5 = fid.variables['xmomentum'][:]
+        ymomentum_5 = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        assert num.allclose(stage_5[0,0], 5*first_value/100)  #Meters
+        assert num.allclose(stage_5[0,3], 5*fourth_value/100)  #Meters
+
+        assert num.allclose(5*stage_1, stage_5)
+
+        # Momentum will also be changed due to new depth
+
+        depth_1 = stage_1-elevation
+        depth_5 = stage_5-elevation
+
+
+        for i in range(stage_1.shape[0]):
+            for j in range(stage_1.shape[1]):            
+                if depth_1[i,j] > epsilon:
+
+                    scale = depth_5[i,j]/depth_1[i,j]
+                    ref_xmomentum = xmomentum_1[i,j] * scale
+                    ref_ymomentum = ymomentum_1[i,j] * scale
+                    
+                    #print i, scale, xmomentum_1[i,j], xmomentum_5[i,j]
+                    
+                    assert num.allclose(xmomentum_5[i,j], ref_xmomentum)
+                    assert num.allclose(ymomentum_5[i,j], ref_ymomentum)
+                    
+        
+
+        fid.close()
+
+
+        #Cleanup
+        import os
+        os.remove(self.test_MOST_file + '.sww')
+
+
+
+    def test_ferret2sww_2(self):
+        """Test that georeferencing etc works when converting from
+        ferret format (lat/lon) to sww format (UTM)
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
+        #
+        # First value (index=0) in small_ha.nc is 0.3400644 cm,
+        # Fourth value (index==3) is -6.50198 cm
+
+
+        from anuga.coordinate_transforms.redfearn import redfearn
+        #fid = NetCDFFile('small_ha.nc')
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+
+        #Pick a coordinate and a value
+
+        time_index = 1
+        lat_index = 0
+        lon_index = 2
+
+        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
+        test_time = fid.variables['TIME'][:][time_index]
+        test_lat = fid.variables['LAT'][:][lat_index]
+        test_lon = fid.variables['LON'][:][lon_index]
+
+        linear_point_index = lat_index*4 + lon_index
+        fid.close()
+
+        #Call conversion (with zero origin)
+        ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 0, 0))
+
+
+        #Work out the UTM coordinates for test point
+        zone, e, n = redfearn(test_lat, test_lon)
+
+        #Read output file 'small.sww'
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        #Check that test coordinate is correctly represented
+        assert num.allclose(x[linear_point_index], e)
+        assert num.allclose(y[linear_point_index], n)
+
+        #Check test value
+        stage = fid.variables['stage'][:]
+
+        assert num.allclose(stage[time_index, linear_point_index], test_value/100)
+
+        fid.close()
+
+        #Cleanup
+        import os
+        os.remove(self.test_MOST_file + '.sww')
+
+
+    def test_ferret2sww_lat_long(self):
+        # Test that min lat long works
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        
+        #Read
+        from anuga.coordinate_transforms.redfearn import redfearn
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+        first_value = fid.variables['HA'][:][0,0,0]
+        fourth_value = fid.variables['HA'][:][0,0,3]
+        fid.close()
+
+
+        #Call conversion (with zero origin)
+        #ferret2sww('small', verbose=self.verbose,
+        #           origin = (56, 0, 0))
+        ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 0, 0), minlat=-34.5, maxlat=-34)
+
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)
+        #print zone, e, n
+
+        #Read output file 'small.sww'
+        #fid = NetCDFFile('small.sww')
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        #Check that first coordinate is correctly represented
+        assert 16 == len(x)
+
+        fid.close()
+
+        #Cleanup
+        import os
+        os.remove(self.test_MOST_file + '.sww')
+
+
+    def test_ferret2sww_lat_longII(self):
+        # Test that min lat long works
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        
+        #Read
+        from anuga.coordinate_transforms.redfearn import redfearn
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+        first_value = fid.variables['HA'][:][0,0,0]
+        fourth_value = fid.variables['HA'][:][0,0,3]
+        fid.close()
+
+
+        #Call conversion (with zero origin)
+        #ferret2sww('small', verbose=False,
+        #           origin = (56, 0, 0))
+        ferret2sww(self.test_MOST_file, verbose=False,
+                   origin = (56, 0, 0), minlat=-34.4, maxlat=-34.2)
+
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)
+        #print zone, e, n
+
+        #Read output file 'small.sww'
+        #fid = NetCDFFile('small.sww')
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        #Check that first coordinate is correctly represented
+        assert 12 == len(x)
+
+        fid.close()
+
+        #Cleanup
+        import os
+        os.remove(self.test_MOST_file + '.sww')
+        
+    def test_ferret2sww3(self):
+        """Elevation included
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        # ELEVATION = [-1 -2 -3 -4
+        #             -5 -6 -7 -8
+        #              ...
+        #              ...      -16]
+        # where the top left corner is -1m,
+        # and the ll corner is -13.0m
+        #
+        # First value (index=0) in small_ha.nc is 0.3400644 cm,
+        # Fourth value (index==3) is -6.50198 cm
+
+        from anuga.coordinate_transforms.redfearn import redfearn
+        import os
+        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
+        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
+        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
+        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
+
+        h1_list = [150.66667,150.83334,151.]
+        h2_list = [-34.5,-34.33333]
+
+        long_name = 'LON'
+        lat_name = 'LAT'
+        time_name = 'TIME'
+
+        nx = 3
+        ny = 2
+
+        for fid in [fid1,fid2,fid3]:
+            fid.createDimension(long_name,nx)
+            fid.createVariable(long_name,netcdf_float,(long_name,))
+            fid.variables[long_name].point_spacing='uneven'
+            fid.variables[long_name].units='degrees_east'
+            fid.variables[long_name].assignValue(h1_list)
+
+            fid.createDimension(lat_name,ny)
+            fid.createVariable(lat_name,netcdf_float,(lat_name,))
+            fid.variables[lat_name].point_spacing='uneven'
+            fid.variables[lat_name].units='degrees_north'
+            fid.variables[lat_name].assignValue(h2_list)
+
+            fid.createDimension(time_name,2)
+            fid.createVariable(time_name,netcdf_float,(time_name,))
+            fid.variables[time_name].point_spacing='uneven'
+            fid.variables[time_name].units='seconds'
+            fid.variables[time_name].assignValue([0.,1.])
+            #if fid == fid3: break
+
+
+        for fid in [fid4]:
+            fid.createDimension(long_name,nx)
+            fid.createVariable(long_name,netcdf_float,(long_name,))
+            fid.variables[long_name].point_spacing='uneven'
+            fid.variables[long_name].units='degrees_east'
+            fid.variables[long_name].assignValue(h1_list)
+
+            fid.createDimension(lat_name,ny)
+            fid.createVariable(lat_name,netcdf_float,(lat_name,))
+            fid.variables[lat_name].point_spacing='uneven'
+            fid.variables[lat_name].units='degrees_north'
+            fid.variables[lat_name].assignValue(h2_list)
+
+        name = {}
+        name[fid1]='HA'
+        name[fid2]='UA'
+        name[fid3]='VA'
+        name[fid4]='ELEVATION'
+
+        units = {}
+        units[fid1]='cm'
+        units[fid2]='cm/s'
+        units[fid3]='cm/s'
+        units[fid4]='m'
+
+        values = {}
+        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
+        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
+        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
+        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
+
+        for fid in [fid1,fid2,fid3]:
+          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
+          fid.variables[name[fid]].point_spacing='uneven'
+          fid.variables[name[fid]].units=units[fid]
+          fid.variables[name[fid]].assignValue(values[fid])
+          fid.variables[name[fid]].missing_value = -99999999.
+          #if fid == fid3: break
+
+        for fid in [fid4]:
+            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
+            fid.variables[name[fid]].point_spacing='uneven'
+            fid.variables[name[fid]].units=units[fid]
+            fid.variables[name[fid]].assignValue(values[fid])
+            fid.variables[name[fid]].missing_value = -99999999.
+
+
+        fid1.sync(); fid1.close()
+        fid2.sync(); fid2.close()
+        fid3.sync(); fid3.close()
+        fid4.sync(); fid4.close()
+
+        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
+        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
+        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
+
+
+        first_amp = fid1.variables['HA'][:][0,0,0]
+        third_amp = fid1.variables['HA'][:][0,0,2]
+        first_elevation = fid2.variables['ELEVATION'][0,0]
+        third_elevation= fid2.variables['ELEVATION'][:][0,2]
+        first_speed = fid3.variables['VA'][0,0,0]
+        third_speed = fid3.variables['VA'][:][0,0,2]
+
+        fid1.close()
+        fid2.close()
+        fid3.close()
+
+        #Call conversion (with zero origin)
+        ferret2sww('test', verbose=self.verbose,
+                   origin = (56, 0, 0), inverted_bathymetry=False)
+
+        os.remove('test_va.nc')
+        os.remove('test_ua.nc')
+        os.remove('test_ha.nc')
+        os.remove('test_e.nc')
+
+        #Read output file 'test.sww'
+        fid = NetCDFFile('test.sww')
+
+
+        #Check first value
+        elevation = fid.variables['elevation'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
+        first_height = first_amp/100 - first_elevation
+        third_height = third_amp/100 - third_elevation
+        first_momentum=first_speed*first_height/100
+        third_momentum=third_speed*third_height/100
+
+        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
+        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
+
+        fid.close()
+
+        #Cleanup
+        os.remove('test.sww')
+
+
+
+    def test_ferret2sww4(self):
+        """Like previous but with augmented variable names as
+        in files produced by ferret as opposed to MOST
+        """
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        # ELEVATION = [-1 -2 -3 -4
+        #             -5 -6 -7 -8
+        #              ...
+        #              ...      -16]
+        # where the top left corner is -1m,
+        # and the ll corner is -13.0m
+        #
+        # First value (index=0) in small_ha.nc is 0.3400644 cm,
+        # Fourth value (index==3) is -6.50198 cm
+
+        from anuga.coordinate_transforms.redfearn import redfearn
+        import os
+        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_w)
+        fid2 = NetCDFFile('test_ua.nc',netcdf_mode_w)
+        fid3 = NetCDFFile('test_va.nc',netcdf_mode_w)
+        fid4 = NetCDFFile('test_e.nc',netcdf_mode_w)
+
+        h1_list = [150.66667,150.83334,151.]
+        h2_list = [-34.5,-34.33333]
+
+#        long_name = 'LON961_1261'
+#        lat_name = 'LAT481_841'
+#        time_name = 'TIME1'
+
+        long_name = 'LON'
+        lat_name = 'LAT'
+        time_name = 'TIME'
+
+        nx = 3
+        ny = 2
+
+        for fid in [fid1,fid2,fid3]:
+            fid.createDimension(long_name,nx)
+            fid.createVariable(long_name,netcdf_float,(long_name,))
+            fid.variables[long_name].point_spacing='uneven'
+            fid.variables[long_name].units='degrees_east'
+            fid.variables[long_name].assignValue(h1_list)
+
+            fid.createDimension(lat_name,ny)
+            fid.createVariable(lat_name,netcdf_float,(lat_name,))
+            fid.variables[lat_name].point_spacing='uneven'
+            fid.variables[lat_name].units='degrees_north'
+            fid.variables[lat_name].assignValue(h2_list)
+
+            fid.createDimension(time_name,2)
+            fid.createVariable(time_name,netcdf_float,(time_name,))
+            fid.variables[time_name].point_spacing='uneven'
+            fid.variables[time_name].units='seconds'
+            fid.variables[time_name].assignValue([0.,1.])
+            #if fid == fid3: break
+
+
+        for fid in [fid4]:
+            fid.createDimension(long_name,nx)
+            fid.createVariable(long_name,netcdf_float,(long_name,))
+            fid.variables[long_name].point_spacing='uneven'
+            fid.variables[long_name].units='degrees_east'
+            fid.variables[long_name].assignValue(h1_list)
+
+            fid.createDimension(lat_name,ny)
+            fid.createVariable(lat_name,netcdf_float,(lat_name,))
+            fid.variables[lat_name].point_spacing='uneven'
+            fid.variables[lat_name].units='degrees_north'
+            fid.variables[lat_name].assignValue(h2_list)
+
+        name = {}
+        name[fid1]='HA'
+        name[fid2]='UA'
+        name[fid3]='VA'
+        name[fid4]='ELEVATION'
+
+        units = {}
+        units[fid1]='cm'
+        units[fid2]='cm/s'
+        units[fid3]='cm/s'
+        units[fid4]='m'
+
+        values = {}
+        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
+        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
+        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
+        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
+
+        for fid in [fid1,fid2,fid3]:
+          fid.createVariable(name[fid],netcdf_float,(time_name,lat_name,long_name))
+          fid.variables[name[fid]].point_spacing='uneven'
+          fid.variables[name[fid]].units=units[fid]
+          fid.variables[name[fid]].assignValue(values[fid])
+          fid.variables[name[fid]].missing_value = -99999999.
+          #if fid == fid3: break
+
+        for fid in [fid4]:
+            fid.createVariable(name[fid],netcdf_float,(lat_name,long_name))
+            fid.variables[name[fid]].point_spacing='uneven'
+            fid.variables[name[fid]].units=units[fid]
+            fid.variables[name[fid]].assignValue(values[fid])
+            fid.variables[name[fid]].missing_value = -99999999.
+
+
+        fid1.sync(); fid1.close()
+        fid2.sync(); fid2.close()
+        fid3.sync(); fid3.close()
+        fid4.sync(); fid4.close()
+
+        fid1 = NetCDFFile('test_ha.nc',netcdf_mode_r)
+        fid2 = NetCDFFile('test_e.nc',netcdf_mode_r)
+        fid3 = NetCDFFile('test_va.nc',netcdf_mode_r)
+
+
+        first_amp = fid1.variables['HA'][:][0,0,0]
+        third_amp = fid1.variables['HA'][:][0,0,2]
+        first_elevation = fid2.variables['ELEVATION'][0,0]
+        third_elevation= fid2.variables['ELEVATION'][:][0,2]
+        first_speed = fid3.variables['VA'][0,0,0]
+        third_speed = fid3.variables['VA'][:][0,0,2]
+
+        fid1.close()
+        fid2.close()
+        fid3.close()
+
+        #Call conversion (with zero origin)
+        ferret2sww('test', verbose=self.verbose, origin = (56, 0, 0)
+                   , inverted_bathymetry=False)
+
+        os.remove('test_va.nc')
+        os.remove('test_ua.nc')
+        os.remove('test_ha.nc')
+        os.remove('test_e.nc')
+
+        #Read output file 'test.sww'
+        fid = NetCDFFile('test.sww')
+
+
+        #Check first value
+        elevation = fid.variables['elevation'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
+        first_height = first_amp/100 - first_elevation
+        third_height = third_amp/100 - third_elevation
+        first_momentum=first_speed*first_height/100
+        third_momentum=third_speed*third_height/100
+
+        assert num.allclose(ymomentum[0][0],first_momentum)  #Meters
+        assert num.allclose(ymomentum[0][2],third_momentum)  #Meters
+
+        fid.close()
+
+        #Cleanup
+        os.remove('test.sww')
+
+
+
+
+    def test_ferret2sww_nz_origin(self):
+        from Scientific.IO.NetCDF import NetCDFFile
+        from anuga.coordinate_transforms.redfearn import redfearn
+
+        #Call conversion (with nonzero origin)
+        ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 100000, 200000))
+
+
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)
+
+        #Read output file 'small.sww'
+        #fid = NetCDFFile('small.sww', netcdf_mode_r)
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        #Check that first coordinate is correctly represented
+        assert num.allclose(x[0], e-100000)
+        assert num.allclose(y[0], n-200000)
+
+        fid.close()
+
+        #Cleanup
+        os.remove(self.test_MOST_file + '.sww')
+
+
+    def test_ferret2sww_lat_longII(self):
+        # Test that min lat long works
+
+        #The test file has
+        # LON = 150.66667, 150.83334, 151, 151.16667
+        # LAT = -34.5, -34.33333, -34.16667, -34 ;
+        
+        #Read
+        from anuga.coordinate_transforms.redfearn import redfearn
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
+        first_value = fid.variables['HA'][:][0,0,0]
+        fourth_value = fid.variables['HA'][:][0,0,3]
+        fid.close()
+
+
+        #Call conversion (with zero origin)
+        #ferret2sww('small', verbose=self.verbose,
+        #           origin = (56, 0, 0))
+        try:
+            ferret2sww(self.test_MOST_file, verbose=self.verbose,
+                   origin = (56, 0, 0), minlat=-34.5, maxlat=-35)
+        except AssertionError:
+            pass
+        else:
+            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
+
+    def test_sww_extent(self):
+        """Not a test, rather a look at the sww format
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        self.domain.set_name('datatest' + str(id(self)))
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.reduction = mean
+        self.domain.set_datadir('.')
+        #self.domain.tight_slope_limiters = 1        
+
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.time = 2.
+
+        #Modify stage at second timestep
+        stage = self.domain.quantities['stage'].vertex_values
+        self.domain.set_quantity('stage', stage/2)
+
+        sww.store_timestep()
+
+        file_and_extension_name = self.domain.get_name() + ".sww"
+        #print "file_and_extension_name",file_and_extension_name
+        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
+               extent_sww(file_and_extension_name )
+
+        assert num.allclose(xmin, 0.0)
+        assert num.allclose(xmax, 1.0)
+        assert num.allclose(ymin, 0.0)
+        assert num.allclose(ymax, 1.0)
+
+        # FIXME (Ole): Revisit these numbers
+        #assert num.allclose(stagemin, -0.85), 'stagemin=%.4f' %stagemin
+        #assert num.allclose(stagemax, 0.15), 'stagemax=%.4f' %stagemax
+
+
+        #Cleanup
+        os.remove(sww.filename)
+
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_File_Conversion, 'test_sww')
+    
+    # FIXME(Ole): When Ross has implemented logging, we can 
+    # probably get rid of all this:
+    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
+        Test_File_Conversion.verbose=True
+        saveout = sys.stdout   
+        filename = ".temp_verbose"
+        fid = open(filename, 'w')
+        sys.stdout = fid
+    else:
+        pass
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
+
+    # Cleaning up
+    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
+        sys.stdout = saveout 
+        fid.close() 
+        os.remove(filename)
+
+
+    def test_asc_csiro2sww(self):
+        import tempfile
+
+        bath_dir = tempfile.mkdtemp()
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
+        #bath_dir = 'bath_data_manager_test'
+        #print "os.getcwd( )",os.getcwd( )
+        elevation_dir =  tempfile.mkdtemp()
+        #elevation_dir = 'elev_expanded'
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
+        fid = open(bath_dir_filename, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    0.000    3000.0
+     0.000     9000.000     0.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000
+""")
+        fid.close()
+
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+
+        fid = open(ucur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(ucur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+
+        fid = open(vcur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(vcur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        sww_file = 'a_test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file,
+                      verbose=self.verbose)
+
+        # check the sww file
+
+        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        geo_ref = Geo_reference(NetCDFObject=fid)
+        #print "geo_ref",geo_ref
+        x_ref = geo_ref.get_xllcorner()
+        y_ref = geo_ref.get_yllcorner()
+        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
+        assert num.allclose(x_ref, 587798.418) # (-38, 148)
+        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
+
+        #Zone:   55
+        #Easting:  588095.674  Northing: 5821451.722
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
+        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
+
+        #Zone:   55
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
+        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
+
+        #Zone:   55
+        #Easting:  609748.788  Northing: 5793447.860
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
+        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
+
+        assert num.allclose(z[0],9000.0 )
+        assert num.allclose(stage[0][1],0.0 )
+
+        #(4000+1000)*60
+        assert num.allclose(xmomentum[1][1],300000.0 )
+
+
+        fid.close()
+
+        #tidy up
+        os.remove(bath_dir_filename)
+        os.rmdir(bath_dir)
+
+        os.remove(elevation_dir_filename1)
+        os.remove(elevation_dir_filename2)
+        os.rmdir(elevation_dir)
+
+        os.remove(ucur_dir_filename1)
+        os.remove(ucur_dir_filename2)
+        os.rmdir(ucur_dir)
+
+        os.remove(vcur_dir_filename1)
+        os.remove(vcur_dir_filename2)
+        os.rmdir(vcur_dir)
+
+
+        # remove sww file
+        os.remove(sww_file)
+
+    def test_asc_csiro2sww2(self):
+        import tempfile
+
+        bath_dir = tempfile.mkdtemp()
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
+        #bath_dir = 'bath_data_manager_test'
+        #print "os.getcwd( )",os.getcwd( )
+        elevation_dir =  tempfile.mkdtemp()
+        #elevation_dir = 'elev_expanded'
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
+        fid = open(bath_dir_filename, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    0.000    3000.0
+     0.000     -9999.000     -9999.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000
+""")
+        fid.close()
+
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+
+        fid = open(ucur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(ucur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+
+        fid = open(vcur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(vcur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        try:
+            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
+                          vcur_dir, sww_file,
+                      verbose=self.verbose)
+        except:
+            #tidy up
+            os.remove(bath_dir_filename)
+            os.rmdir(bath_dir)
+
+            os.remove(elevation_dir_filename1)
+            os.remove(elevation_dir_filename2)
+            os.rmdir(elevation_dir)
+
+            os.remove(ucur_dir_filename1)
+            os.remove(ucur_dir_filename2)
+            os.rmdir(ucur_dir)
+
+            os.remove(vcur_dir_filename1)
+            os.remove(vcur_dir_filename2)
+            os.rmdir(vcur_dir)
+        else:
+            #tidy up
+            os.remove(bath_dir_filename)
+            os.rmdir(bath_dir)
+
+            os.remove(elevation_dir_filename1)
+            os.remove(elevation_dir_filename2)
+            os.rmdir(elevation_dir)
+            raise 'Should raise exception'
+
+            os.remove(ucur_dir_filename1)
+            os.remove(ucur_dir_filename2)
+            os.rmdir(ucur_dir)
+
+            os.remove(vcur_dir_filename1)
+            os.remove(vcur_dir_filename2)
+            os.rmdir(vcur_dir)
+
+
+
+    def test_asc_csiro2sww3(self):
+        import tempfile
+
+        bath_dir = tempfile.mkdtemp()
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
+        #bath_dir = 'bath_data_manager_test'
+        #print "os.getcwd( )",os.getcwd( )
+        elevation_dir =  tempfile.mkdtemp()
+        #elevation_dir = 'elev_expanded'
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
+        fid = open(bath_dir_filename, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    0.000    3000.0
+     0.000     -9999.000     -9999.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000
+""")
+        fid.close()
+
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+
+        fid = open(ucur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(ucur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+
+        fid = open(vcur_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+        fid = open(vcur_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000
+""")
+        fid.close()
+
+        sww_file = 'a_test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
+                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
+                      mean_stage = 100,
+                      verbose=self.verbose)
+
+        # check the sww file
+
+        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        geo_ref = Geo_reference(NetCDFObject=fid)
+        #print "geo_ref",geo_ref
+        x_ref = geo_ref.get_xllcorner()
+        y_ref = geo_ref.get_yllcorner()
+        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
+        assert num.allclose(x_ref, 587798.418) # (-38, 148)
+        assert num.allclose(y_ref, 5793123.477)# (-38, 148.5)
+
+        #Zone:   55
+        #Easting:  588095.674  Northing: 5821451.722
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
+        assert num.allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
+
+        #Zone:   55
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
+        assert num.allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
+
+        #Zone:   55
+        #Easting:  609748.788  Northing: 5793447.860
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
+        assert num.allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
+
+        assert num.allclose(z[0],9000.0 )
+        assert num.allclose(stage[0][4],100.0 )
+        assert num.allclose(stage[0][5],100.0 )
+
+        #(100.0 - 9000)*10
+        assert num.allclose(xmomentum[0][4], -89000.0 )
+
+        #(100.0 - -1000.000)*10
+        assert num.allclose(xmomentum[0][5], 11000.0 )
+
+        fid.close()
+
+        #tidy up
+        os.remove(bath_dir_filename)
+        os.rmdir(bath_dir)
+
+        os.remove(elevation_dir_filename1)
+        os.remove(elevation_dir_filename2)
+        os.rmdir(elevation_dir)
+
+        os.remove(ucur_dir_filename1)
+        os.remove(ucur_dir_filename2)
+        os.rmdir(ucur_dir)
+
+        os.remove(vcur_dir_filename1)
+        os.remove(vcur_dir_filename2)
+        os.rmdir(vcur_dir)
+
+        # remove sww file
+        os.remove(sww_file)
+
+
+    def test_asc_csiro2sww4(self):
+        """
+        Test specifying the extent
+        """
+
+        import tempfile
+
+        bath_dir = tempfile.mkdtemp()
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
+        #bath_dir = 'bath_data_manager_test'
+        #print "os.getcwd( )",os.getcwd( )
+        elevation_dir =  tempfile.mkdtemp()
+        #elevation_dir = 'elev_expanded'
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
+        fid = open(bath_dir_filename, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -9000.000    -1000.000   -3000.0 -2000.000
+   -1000.000    9000.000  -1000.000 -3000.000
+   -4000.000    6000.000   2000.000 -5000.000
+   -9000.000    -1000.000   -3000.0 -2000.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename1, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -900.000    -100.000   -300.0 -200.000
+   -100.000    900.000  -100.000 -300.000
+   -400.000    600.000   200.000 -500.000
+   -900.000    -100.000   -300.0 -200.000
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename2, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -990.000    -110.000   -330.0 -220.000
+   -110.000    990.000  -110.000 -330.000
+   -440.000    660.000   220.000 -550.000
+   -990.000    -110.000   -330.0 -220.000
+""")
+        fid.close()
+
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+
+        fid = open(ucur_dir_filename1, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -90.000    -10.000   -30.0 -20.000
+   -10.000    90.000  -10.000 -30.000
+   -40.000    60.000   20.000 -50.000
+   -90.000    -10.000   -30.0 -20.000
+""")
+        fid.close()
+        fid = open(ucur_dir_filename2, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -90.000    -10.000   -30.0 -20.000
+   -10.000    99.000  -11.000 -30.000
+   -40.000    66.000   22.000 -50.000
+   -90.000    -10.000   -30.0 -20.000
+""")
+        fid.close()
+
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+
+        fid = open(vcur_dir_filename1, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -90.000    -10.000   -30.0 -20.000
+   -10.000    80.000  -20.000 -30.000
+   -40.000    50.000   10.000 -50.000
+   -90.000    -10.000   -30.0 -20.000
+""")
+        fid.close()
+        fid = open(vcur_dir_filename2, 'w')
+        fid.write(""" ncols             4
+ nrows             4
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+   -90.000    -10.000   -30.0 -20.000
+   -10.000    88.000  -22.000 -30.000
+   -40.000    55.000   11.000 -50.000
+   -90.000    -10.000   -30.0 -20.000
+""")
+        fid.close()
+
+        sww_file = tempfile.mktemp(".sww")
+        #sww_file = 'a_test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
+                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
+                      mean_stage = 100,
+                       minlat = -37.6, maxlat = -37.6,
+                  minlon = 148.3, maxlon = 148.3,
+                      verbose=self.verbose
+                      #,verbose = True
+                      )
+
+        # check the sww file
+
+        fid = NetCDFFile(sww_file, netcdf_mode_r)    # Open existing file for read
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        geo_ref = Geo_reference(NetCDFObject=fid)
+        #print "geo_ref",geo_ref
+        x_ref = geo_ref.get_xllcorner()
+        y_ref = geo_ref.get_yllcorner()
+        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
+
+        assert num.allclose(fid.starttime, 0.0) # (-37.45, 148.25)
+        assert num.allclose(x_ref, 610120.388) # (-37.45, 148.25)
+        assert num.allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
+
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
+
+        #print "x",x
+        #print "y",y
+        self.failUnless(len(x) == 4,'failed') # 2*2
+        self.failUnless(len(x) == 4,'failed') # 2*2
+
+        #Zone:   55
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
+        # magic number - y is close enough for me.
+        assert num.allclose(x[3], 632145.63 - x_ref)
+        assert num.allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
+
+        assert num.allclose(z[0],9000.0 ) #z is elevation info
+        #print "z",z
+        # 2 time steps, 4 points
+        self.failUnless(xmomentum.shape == (2,4), 'failed')
+        self.failUnless(ymomentum.shape == (2,4), 'failed')
+
+        #(100.0 - -1000.000)*10
+        #assert num.allclose(xmomentum[0][5], 11000.0 )
+
+        fid.close()
+
+        # is the sww file readable?
+        #Lets see if we can convert it to a dem!
+        # if you uncomment, remember to delete the file
+        #print "sww_file",sww_file
+        #dem_file = tempfile.mktemp(".dem")
+        domain = sww2domain(sww_file) ###, dem_file)
+        domain.check_integrity()
+
+        #tidy up
+        os.remove(bath_dir_filename)
+        os.rmdir(bath_dir)
+
+        os.remove(elevation_dir_filename1)
+        os.remove(elevation_dir_filename2)
+        os.rmdir(elevation_dir)
+
+        os.remove(ucur_dir_filename1)
+        os.remove(ucur_dir_filename2)
+        os.rmdir(ucur_dir)
+
+        os.remove(vcur_dir_filename1)
+        os.remove(vcur_dir_filename2)
+        os.rmdir(vcur_dir)
+
+        # remove sww file
+        os.remove(sww_file)
+
+
+
+    def test_get_min_max_indexes(self):
+        latitudes = [3,2,1,0]
+        longitudes = [0,10,20,30]
+
+        # k - lat
+        # l - lon
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            -10,4,-10,31)
+
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+        self.failUnless(latitudes == latitudes_new and \
+                        longitudes == longitudes_news,
+                         'failed')
+
+        ## 2nd test
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            0.5,2.5,5,25)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes == latitudes_new and \
+                        longitudes == longitudes_news,
+                         'failed')
+
+        ## 3rd test
+        kmin, kmax, lmin, lmax = get_min_max_indexes(\
+            latitudes,
+            longitudes,
+            1.1,1.9,12,17)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [2, 1] and \
+                        longitudes_news == [10, 20],
+                         'failed')
+
+
+        ## 4th test
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+                                                      -0.1,1.9,-2,17)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [2, 1, 0] and \
+                        longitudes_news == [0, 10, 20],
+                         'failed')
+        ## 5th test
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            0.1,1.9,2,17)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [2, 1, 0] and \
+                        longitudes_news == [0, 10, 20],
+                         'failed')
+
+        ## 6th test
+
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            1.5,4,18,32)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [3, 2, 1] and \
+                        longitudes_news == [10, 20, 30],
+                         'failed')
+
+
+        ## 7th test
+        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            1.5,1.5,15,15)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        m2d = m2d[kmin:kmax,lmin:lmax]
+        #print "m2d", m2d
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(num.alltrue(latitudes_new == [2, 1]) and 
+                        num.alltrue(longitudes_news == [10, 20]),
+                        'failed')
+
+        self.failUnless(num.alltrue(m2d == [[5,6],[9,10]]), 'failed')
+
+    def test_get_min_max_indexes_lat_ascending(self):
+        latitudes = [0,1,2,3]
+        longitudes = [0,10,20,30]
+
+        # k - lat
+        # l - lon
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            -10,4,-10,31)
+
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+        self.failUnless(latitudes == latitudes_new and \
+                        longitudes == longitudes_news,
+                         'failed')
+
+        ## 3rd test
+        kmin, kmax, lmin, lmax = get_min_max_indexes(\
+            latitudes,
+            longitudes,
+            1.1,1.9,12,17)
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [1, 2] and \
+                        longitudes_news == [10, 20],
+                         'failed')
+
+    def test_get_min_max_indexes2(self):
+        latitudes = [-30,-35,-40,-45]
+        longitudes = [148,149,150,151]
+
+        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.int)    #array default#
+
+        # k - lat
+        # l - lon
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            -37,-27,147,149.5)
+
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        #print "m2d", m2d
+        #print "latitudes", latitudes
+        #print "longitudes",longitudes
+        #print "latitudes[kmax]", latitudes[kmax]
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_new = longitudes[lmin:lmax]
+        m2d = m2d[kmin:kmax,lmin:lmax]
+        #print "m2d", m2d
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_new",longitudes_new
+
+        self.failUnless(latitudes_new == [-30, -35, -40] and
+                        longitudes_new == [148, 149,150],
+                        'failed')
+        self.failUnless(num.alltrue(m2d == [[0,1,2],[4,5,6],[8,9,10]]),
+                        'failed')
+
+    def test_get_min_max_indexes3(self):
+        latitudes = [-30,-35,-40,-45,-50,-55,-60]
+        longitudes = [148,149,150,151]
+
+        # k - lat
+        # l - lon
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes,
+            -43,-37,148.5,149.5)
+
+
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        #print "latitudes", latitudes
+        #print "longitudes",longitudes
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == [-35, -40, -45] and
+                        longitudes_news == [148, 149,150],
+                         'failed')
+
+    def test_get_min_max_indexes4(self):
+        latitudes = [-30,-35,-40,-45,-50,-55,-60]
+        longitudes = [148,149,150,151]
+
+        # k - lat
+        # l - lon
+        kmin, kmax, lmin, lmax = get_min_max_indexes(
+            latitudes,longitudes)
+
+        #print "kmin",kmin;print "kmax",kmax
+        #print "lmin",lmin;print "lmax",lmax
+        #print "latitudes", latitudes
+        #print "longitudes",longitudes
+        latitudes_new = latitudes[kmin:kmax]
+        longitudes_news = longitudes[lmin:lmax]
+        #print "latitudes_new", latitudes_new
+        #print "longitudes_news",longitudes_news
+
+        self.failUnless(latitudes_new == latitudes  and \
+                        longitudes_news == longitudes,
+                         'failed')
+
+    def test_tsh2sww(self):
+        import os
+        import tempfile
+
+        tsh_file = tempfile.mktemp(".tsh")
+        file = open(tsh_file,"w")
+        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
+0 0.0 0.0 0.0 0.0 0.01 \n \
+1 1.0 0.0 10.0 10.0 0.02  \n \
+2 0.0 1.0 0.0 10.0 0.03  \n \
+3 0.5 0.25 8.0 12.0 0.04  \n \
+# Vert att title  \n \
+elevation  \n \
+stage  \n \
+friction  \n \
+2 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
+0 0 3 2 -1  -1  1 dsg\n\
+1 0 1 3 -1  0 -1   ole nielsen\n\
+4 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
+0 1 0 2 \n\
+1 0 2 3 \n\
+2 2 3 \n\
+3 3 1 1 \n\
+3 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
+0 216.0 -86.0 \n \
+1 160.0 -167.0 \n \
+2 114.0 -91.0 \n \
+3 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
+0 0 1 0 \n \
+1 1 2 0 \n \
+2 2 0 0 \n \
+0 # <x> <y> ...Mesh Holes... \n \
+0 # <x> <y> <attribute>...Mesh Regions... \n \
+0 # <x> <y> <attribute>...Mesh Regions, area... \n\
+#Geo reference \n \
+56 \n \
+140 \n \
+120 \n")
+        file.close()
+
+        #sww_file = tempfile.mktemp(".sww")
+        #print "sww_file",sww_file
+        #print "sww_file",tsh_file
+        tsh2sww(tsh_file,
+                verbose=self.verbose)
+
+        os.remove(tsh_file)
+        os.remove(tsh_file[:-4] + '.sww')
+        
+          
+    def test_urs2sww_test_fail(self):
+        points_num = -100
+        time_step_count = 45
+        time_step = -7
+        file_handle, base_name = tempfile.mkstemp("")        
+        os.close(file_handle)
+        os.remove(base_name)
+        
+        files = []
+        quantities = ['HA','UA','VA']
+        
+        mux_names = [WAVEHEIGHT_MUX_LABEL,
+                     EAST_VELOCITY_LABEL,
+                     NORTH_VELOCITY_LABEL]
+        for i,q in enumerate(quantities): 
+            #Write C files
+            columns = 3 # long, lat , depth
+            file = base_name + mux_names[i]
+            f = open(file, 'wb')
+            files.append(file)
+            f.write(pack('i',points_num))
+            f.write(pack('i',time_step_count))
+            f.write(pack('f',time_step))
+
+            f.close()   
+        tide = 1
+        try:
+            urs2sww(base_name, remove_nc_files=True, mean_stage=tide,
+                      verbose=self.verbose)        
+        except ANUGAError:
+            pass
+        else:
+            self.delete_mux(files)
+            msg = 'Should have raised exception'
+            raise msg
+        sww_file = base_name + '.sww'
+        self.delete_mux(files)
+        
+    def test_urs2sww_test_fail2(self):
+        base_name = 'Harry-high-pants'
+        try:
+            urs2sww(base_name)        
+        except IOError:
+            pass
+        else:
+            self.delete_mux(files)
+            msg = 'Should have raised exception'
+            raise msg
+           
+    def test_urs2sww(self):
+        tide = 1
+        base_name, files = self.create_mux()
+        urs2sww(base_name
+                #, origin=(0,0,0)
+                , mean_stage=tide
+                , remove_nc_files=True,
+                      verbose=self.verbose
+                )
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)
+       
+        assert num.allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
+
+        # Make x and y absolute
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+        # = n*(e+tide+n) based on how I'm writing these files
+        # 
+        answer_x = n*(e+tide+n)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert num.allclose(answer_x, actual_x)  #Meters
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        #momentum = velocity_va *(stage+elevation)
+        # = e*(e+tide+n) based on how I'm writing these files
+        answer_y = e*(e+tide+n) * -1 # talking into account mux file format
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert num.allclose(answer_y, actual_y)  #Meters
+        
+        assert num.allclose(answer_x, actual_x)  #Meters
+        
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs2sww_momentum(self):
+        tide = 1
+        time_step_count = 3
+        time_step = 2
+        #lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,114.5), (-21.,115.)]
+        # This is gridded
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        depth=20
+        ha=2
+        ua=5
+        va=-10 #-ve added to take into account mux file format where south
+               # is positive.
+        base_name, files = self.write_mux(lat_long_points,
+                                          time_step_count, time_step,
+                                          depth=depth,
+                                          ha=ha,
+                                          ua=ua,
+                                          va=va)
+        # write_mux(self,lat_long_points, time_step_count, time_step,
+        #          depth=None, ha=None, ua=None, va=None
+        urs2sww(base_name
+                #, origin=(0,0,0)
+                , mean_stage=tide
+                , remove_nc_files=True,
+                      verbose=self.verbose
+                )
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        #assert allclose(stage[0,0], e + tide)  #Meters
+        #print "xmomentum", xmomentum
+        #print "ymomentum", ymomentum
+        #Check the momentums - ua
+        #momentum = velocity*water height
+        #water height = mux_depth + mux_height +tide
+        #water height = mux_depth + mux_height +tide
+        #momentum = velocity*(mux_depth + mux_height +tide)
+        #
+        
+        answer = 115
+        actual = xmomentum[0,0]
+        assert num.allclose(answer, actual)  #Meters^2/ sec
+        answer = 230
+        actual = ymomentum[0,0]
+        #print "answer",answer
+        #print "actual",actual 
+        assert num.allclose(answer, actual)  #Meters^2/ sec
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+
+        #assert allclose(stage[0], x +tide)  #Meters
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        #assert allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+  
+    def test_urs2sww_origin(self):
+        tide = 1
+        base_name, files = self.create_mux()
+        urs2sww(base_name
+                , origin=(0,0,0)
+                , mean_stage=tide
+                , remove_nc_files=True,
+                      verbose=self.verbose
+                )
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+
+        #  x and y are absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+
+        
+        time = fid.variables['time'][:]
+        #print "time", time
+        assert num.allclose([0.,0.5,1.], time)
+        assert fid.starttime == 0.0
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)       
+       
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert num.allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+ 
+    def test_urs2sww_minmaxlatlong(self):
+        
+        #longitudes = [150.66667, 150.83334, 151., 151.16667]
+        #latitudes = [-34.5, -34.33333, -34.16667, -34]
+
+        tide = 1
+        base_name, files = self.create_mux()
+        urs2sww(base_name,
+                minlat=-34.5,
+                maxlat=-34,
+                minlon= 150.66667,
+                maxlon= 151.16667,
+                mean_stage=tide,
+                remove_nc_files=True,
+                      verbose=self.verbose
+                )
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667) 
+        assert num.allclose([x[0],y[0]], [e,n])
+
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert num.allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert num.allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert num.allclose(stage[0], x +tide)  #Meters
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert num.allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs2sww_minmaxmintmaxt(self):
+        
+        #longitudes = [150.66667, 150.83334, 151., 151.16667]
+        #latitudes = [-34.5, -34.33333, -34.16667, -34]
+
+        tide = 1
+        base_name, files = self.create_mux()
+        
+        urs2sww(base_name,
+                mint=0.25,
+                maxt=0.75,
+                mean_stage=tide,
+                remove_nc_files=True,
+                verbose=self.verbose)
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+
+        
+        time = fid.variables['time'][:]
+        assert num.allclose(time, [0.0]) # the time is relative
+        assert fid.starttime == 0.5
+        
+        fid.close()
+        self.delete_mux(files)
+        #print "sww_file", sww_file
+        os.remove(sww_file)
+  
+
+
+#-------------------------------------------------------------
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_File_Conversion,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(suite)
+
Index: /anuga_core/source/anuga/shallow_water/data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7741)
+++ /anuga_core/source/anuga/shallow_water/data_manager.py	(revision 7742)
@@ -1656,215 +1656,4 @@
     fid.close()
 
-
-##
-# @brief Convert ASC file to DEM file.
-# @param basename_in Stem of input filename.
-# @param basename_out Stem of output filename.
-# @param use_cache ??
-# @param verbose True if this function is to be verbose.
-# @return 
-# @note A PRJ file with same stem basename must exist and is used to fix the
-#       UTM zone, datum, false northings and eastings.
-def convert_dem_from_ascii2netcdf(basename_in, basename_out=None,
-                                  use_cache=False,
-                                  verbose=False):
-    """Read Digital Elevation model from the following ASCII format (.asc)
-
-    Example:
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    Convert basename_in + '.asc' to NetCDF format (.dem)
-    mimicking the ASCII format closely.
-
-    An accompanying file with same basename_in but extension .prj must exist
-    and is used to fix the UTM zone, datum, false northings and eastings.
-
-    The prj format is assumed to be as
-
-    Projection    UTM
-    Zone          56
-    Datum         WGS84
-    Zunits        NO
-    Units         METERS
-    Spheroid      WGS84
-    Xshift        0.0000000000
-    Yshift        10000000.0000000000
-    Parameters
-    """
-
-    kwargs = {'basename_out': basename_out, 'verbose': verbose}
-
-    if use_cache is True:
-        from caching import cache
-        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
-                       dependencies=[basename_in + '.asc',
-                                     basename_in + '.prj'],
-                       verbose=verbose)
-
-    else:
-        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
-
-    return result
-
-
-##
-# @brief Convert an ASC file to a DEM file.
-# @param basename_in Stem of input filename.
-# @param basename_out Stem of output filename.
-# @param verbose True if this function is to be verbose.
-def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
-                                   verbose = False):
-    """Read Digital Elevation model from the following ASCII format (.asc)
-
-    Internal function. See public function convert_dem_from_ascii2netcdf
-    for details.
-    """
-
-    import os
-    from Scientific.IO.NetCDF import NetCDFFile
-
-    root = basename_in
-
-    # Read Meta data
-    if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
-
-    metadatafile = open(root + '.prj')
-    metalines = metadatafile.readlines()
-    metadatafile.close()
-
-    L = metalines[0].strip().split()
-    assert L[0].strip().lower() == 'projection'
-    projection = L[1].strip()                   #TEXT
-
-    L = metalines[1].strip().split()
-    assert L[0].strip().lower() == 'zone'
-    zone = int(L[1].strip())
-
-    L = metalines[2].strip().split()
-    assert L[0].strip().lower() == 'datum'
-    datum = L[1].strip()                        #TEXT
-
-    L = metalines[3].strip().split()
-    assert L[0].strip().lower() == 'zunits'     #IGNORE
-    zunits = L[1].strip()                       #TEXT
-
-    L = metalines[4].strip().split()
-    assert L[0].strip().lower() == 'units'
-    units = L[1].strip()                        #TEXT
-
-    L = metalines[5].strip().split()
-    assert L[0].strip().lower() == 'spheroid'   #IGNORE
-    spheroid = L[1].strip()                     #TEXT
-
-    L = metalines[6].strip().split()
-    assert L[0].strip().lower() == 'xshift'
-    false_easting = float(L[1].strip())
-
-    L = metalines[7].strip().split()
-    assert L[0].strip().lower() == 'yshift'
-    false_northing = float(L[1].strip())
-
-    #Read DEM data
-    datafile = open(basename_in + '.asc')
-
-    if verbose: log.critical('Reading DEM from %s' % (basename_in + '.asc'))
-
-    lines = datafile.readlines()
-    datafile.close()
-
-    if verbose: log.critical('Got %d lines' % len(lines))
-
-    ncols = int(lines[0].split()[1].strip())
-    nrows = int(lines[1].split()[1].strip())
-
-    # Do cellsize (line 4) before line 2 and 3
-    cellsize = float(lines[4].split()[1].strip())
-
-    # Checks suggested by Joaquim Luis
-    # Our internal representation of xllcorner
-    # and yllcorner is non-standard.
-    xref = lines[2].split()
-    if xref[0].strip() == 'xllcorner':
-        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
-    elif xref[0].strip() == 'xllcenter':
-        xllcorner = float(xref[1].strip())
-    else:
-        msg = 'Unknown keyword: %s' % xref[0].strip()
-        raise Exception, msg
-
-    yref = lines[3].split()
-    if yref[0].strip() == 'yllcorner':
-        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
-    elif yref[0].strip() == 'yllcenter':
-        yllcorner = float(yref[1].strip())
-    else:
-        msg = 'Unknown keyword: %s' % yref[0].strip()
-        raise Exception, msg
-
-    NODATA_value = int(lines[5].split()[1].strip())
-
-    assert len(lines) == nrows + 6
-
-    if basename_out == None:
-        netcdfname = root + '.dem'
-    else:
-        netcdfname = basename_out + '.dem'
-
-    if verbose: log.critical('Store to NetCDF file %s' % netcdfname)
-
-    # NetCDF file definition
-    fid = NetCDFFile(netcdfname, netcdf_mode_w)
-
-    #Create new file
-    fid.institution = 'Geoscience Australia'
-    fid.description = 'NetCDF DEM format for compact and portable storage ' \
-                      'of spatial point data'
-
-    fid.ncols = ncols
-    fid.nrows = nrows
-    fid.xllcorner = xllcorner
-    fid.yllcorner = yllcorner
-    fid.cellsize = cellsize
-    fid.NODATA_value = NODATA_value
-
-    fid.zone = zone
-    fid.false_easting = false_easting
-    fid.false_northing = false_northing
-    fid.projection = projection
-    fid.datum = datum
-    fid.units = units
-
-    # dimension definitions
-    fid.createDimension('number_of_rows', nrows)
-    fid.createDimension('number_of_columns', ncols)
-
-    # variable definitions
-    fid.createVariable('elevation', netcdf_float, ('number_of_rows',
-                                                   'number_of_columns'))
-
-    # Get handles to the variables
-    elevation = fid.variables['elevation']
-
-    #Store data
-    n = len(lines[6:])
-    for i, line in enumerate(lines[6:]):
-        fields = line.split()
-        if verbose and i % ((n+10)/10) == 0:
-            log.critical('Processing row %d of %d' % (i, nrows))
-           
-        if len(fields) != ncols:
-            msg = 'Wrong number of columns in file "%s" line %d\n' % (basename_in + '.asc', i)
-            msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
-            raise Exception, msg
-
-        elevation[i, :] = num.array([float(x) for x in fields])
-
-    fid.close()
 
 
Index: /anuga_core/source/anuga/shallow_water/sww_file.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/sww_file.py	(revision 7741)
+++ /anuga_core/source/anuga/shallow_water/sww_file.py	(revision 7742)
@@ -7,4 +7,5 @@
 from anuga.config import max_float
 from anuga.utilities.numerical_tools import ensure_numeric
+import anuga.utilities.log as log
 from Scientific.IO.NetCDF import NetCDFFile
 
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7741)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 7742)
@@ -18,7 +18,7 @@
 from anuga.shallow_water.data_manager import *
 from anuga.shallow_water.sww_file import SWW_file
-from anuga.shallow_water.file_conversion import tsh2sww, \
-                        asc_csiro2sww, pmesh_to_domain_instance, \
-                        dem2pts, _get_min_max_indexes
+from anuga.file_conversion.file_conversion import tsh2sww, \
+                        asc_csiro2sww, pmesh_to_domain_instance
+                        
 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
 from anuga.abstract_2d_finite_volumes.util import file_function
@@ -735,380 +735,4 @@
         #Cleanup
         os.remove(sww.filename)
-
-
-
-    def test_dem2pts_bounding_box_v2(self):
-        """Test conversion from dem in ascii format to native NetCDF format
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Write test asc file
-        root = 'demtest'
-
-        filename = root+'.asc'
-        fid = open(filename, 'w')
-        fid.write("""ncols         10
-nrows         10
-xllcorner     2000
-yllcorner     3000
-cellsize      1
-NODATA_value  -9999
-""")
-        #Create linear function
-        ref_points = []
-        ref_elevation = []
-        x0 = 2000
-        y = 3010
-        yvec = range(10)
-        xvec = range(10)
-        z = -1
-        for i in range(10):
-            y = y - 1
-            for j in range(10):
-                x = x0 + xvec[j]
-                z += 1
-                ref_points.append ([x,y])
-                ref_elevation.append(z)
-                fid.write('%f ' %z)
-            fid.write('\n')
-
-        fid.close()
-
-        #print 'sending pts', ref_points
-        #print 'sending elev', ref_elevation
-
-        #Write prj file with metadata
-        metafilename = root+'.prj'
-        fid = open(metafilename, 'w')
-
-
-        fid.write("""Projection UTM
-Zone 56
-Datum WGS84
-Zunits NO
-Units METERS
-Spheroid WGS84
-Xshift 0.0000000000
-Yshift 10000000.0000000000
-Parameters
-""")
-        fid.close()
-
-        #Convert to NetCDF pts
-        convert_dem_from_ascii2netcdf(root)
-        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
-                northing_min=3003.0, northing_max=3006.0,
-                verbose=self.verbose)
-
-        #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
-        assert fid.xllcorner[0] == 2002.0
-        assert fid.yllcorner[0] == 3003.0
-
-        #create new reference points
-        newz = []
-        newz[0:5] = ref_elevation[32:38]
-        newz[6:11] = ref_elevation[42:48]
-        newz[12:17] = ref_elevation[52:58]
-        newz[18:23] = ref_elevation[62:68]
-        ref_elevation = []
-        ref_elevation = newz
-        ref_points = []
-        x0 = 2002
-        y = 3007
-        yvec = range(4)
-        xvec = range(6)
-        for i in range(4):
-            y = y - 1
-            ynew = y - 3003.0
-            for j in range(6):
-                x = x0 + xvec[j]
-                xnew = x - 2002.0
-                ref_points.append ([xnew,ynew]) #Relative point values
-
-        assert num.allclose(points, ref_points)
-
-        assert num.allclose(elevation, ref_elevation)
-
-        #Cleanup
-        fid.close()
-
-
-        os.remove(root + '.pts')
-        os.remove(root + '.dem')
-        os.remove(root + '.asc')
-        os.remove(root + '.prj')
-
-
-    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
-        """Test conversion from dem in ascii format to native NetCDF format
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Write test asc file
-        root = 'demtest'
-
-        filename = root+'.asc'
-        fid = open(filename, 'w')
-        fid.write("""ncols         10
-nrows         10
-xllcorner     2000
-yllcorner     3000
-cellsize      1
-NODATA_value  -9999
-""")
-        #Create linear function
-        ref_points = []
-        ref_elevation = []
-        x0 = 2000
-        y = 3010
-        yvec = range(10)
-        xvec = range(10)
-        #z = range(100)
-        z = num.zeros(100, num.int)     #array default#
-        NODATA_value = -9999
-        count = -1
-        for i in range(10):
-            y = y - 1
-            for j in range(10):
-                x = x0 + xvec[j]
-                ref_points.append ([x,y])
-                count += 1
-                z[count] = (4*i - 3*j)%13
-                if j == 4: z[count] = NODATA_value #column inside clipping region
-                if j == 8: z[count] = NODATA_value #column outside clipping region
-                if i == 9: z[count] = NODATA_value #row outside clipping region
-                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
-                ref_elevation.append( z[count] )
-                fid.write('%f ' %z[count])
-            fid.write('\n')
-
-        fid.close()
-
-        #print 'sending elev', ref_elevation
-
-        #Write prj file with metadata
-        metafilename = root+'.prj'
-        fid = open(metafilename, 'w')
-
-
-        fid.write("""Projection UTM
-Zone 56
-Datum WGS84
-Zunits NO
-Units METERS
-Spheroid WGS84
-Xshift 0.0000000000
-Yshift 10000000.0000000000
-Parameters
-""")
-        fid.close()
-
-        #Convert to NetCDF pts
-        convert_dem_from_ascii2netcdf(root)
-        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
-                northing_min=3003.0, northing_max=3006.0,
-                verbose=self.verbose)
-
-        #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
-        assert fid.xllcorner[0] == 2002.0
-        assert fid.yllcorner[0] == 3003.0
-
-        #create new reference points
-        newz = num.zeros(19, num.int)       #array default#
-        newz[0:2] = ref_elevation[32:34]
-        newz[2:5] = ref_elevation[35:38]
-        newz[5:7] = ref_elevation[42:44]
-        newz[7] = ref_elevation[45]
-        newz[8] = ref_elevation[47]
-        newz[9:11] = ref_elevation[52:54]
-        newz[11:14] = ref_elevation[55:58]
-        newz[14:16] = ref_elevation[62:64]
-        newz[16:19] = ref_elevation[65:68]
-
-
-        ref_elevation = newz
-        ref_points = []
-        new_ref_points = []
-        x0 = 2002
-        y = 3007
-        yvec = range(4)
-        xvec = range(6)
-        for i in range(4):
-            y = y - 1
-            ynew = y - 3003.0
-            for j in range(6):
-                x = x0 + xvec[j]
-                xnew = x - 2002.0
-                if j <> 2 and (i<>1 or j<>4):
-                    ref_points.append([x,y])
-                    new_ref_points.append ([xnew,ynew])
-
-
-        assert num.allclose(points, new_ref_points)
-        assert num.allclose(elevation, ref_elevation)
-
-        #Cleanup
-        fid.close()
-
-
-        os.remove(root + '.pts')
-        os.remove(root + '.dem')
-        os.remove(root + '.asc')
-        os.remove(root + '.prj')
-
-
-    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
-        """Test conversion from dem in ascii format to native NetCDF format
-        Check missing values on clipping boundary
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Write test asc file
-        root = 'demtest'
-
-        filename = root+'.asc'
-        fid = open(filename, 'w')
-        fid.write("""ncols         10
-nrows         10
-xllcorner     2000
-yllcorner     3000
-cellsize      1
-NODATA_value  -9999
-""")
-        #Create linear function
-        ref_points = []
-        ref_elevation = []
-        x0 = 2000
-        y = 3010
-        yvec = range(10)
-        xvec = range(10)
-        #z = range(100)
-        z = num.zeros(100, num.int)     #array default#
-        NODATA_value = -9999
-        count = -1
-        for i in range(10):
-            y = y - 1
-            for j in range(10):
-                x = x0 + xvec[j]
-                ref_points.append ([x,y])
-                count += 1
-                z[count] = (4*i - 3*j)%13
-                if j == 4: z[count] = NODATA_value #column inside clipping region
-                if j == 8: z[count] = NODATA_value #column outside clipping region
-                if i == 6: z[count] = NODATA_value #row on clipping boundary
-                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
-                ref_elevation.append( z[count] )
-                fid.write('%f ' %z[count])
-            fid.write('\n')
-
-        fid.close()
-
-        #print 'sending elev', ref_elevation
-
-        #Write prj file with metadata
-        metafilename = root+'.prj'
-        fid = open(metafilename, 'w')
-
-
-        fid.write("""Projection UTM
-Zone 56
-Datum WGS84
-Zunits NO
-Units METERS
-Spheroid WGS84
-Xshift 0.0000000000
-Yshift 10000000.0000000000
-Parameters
-""")
-        fid.close()
-
-        #Convert to NetCDF pts
-        convert_dem_from_ascii2netcdf(root)
-        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
-                northing_min=3003.0, northing_max=3006.0,
-                verbose=self.verbose)
-
-        #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
-        assert fid.xllcorner[0] == 2002.0
-        assert fid.yllcorner[0] == 3003.0
-
-        #create new reference points
-        newz = num.zeros(14, num.int)       #array default#
-        newz[0:2] = ref_elevation[32:34]
-        newz[2:5] = ref_elevation[35:38]
-        newz[5:7] = ref_elevation[42:44]
-        newz[7] = ref_elevation[45]
-        newz[8] = ref_elevation[47]
-        newz[9:11] = ref_elevation[52:54]
-        newz[11:14] = ref_elevation[55:58]
-
-
-
-        ref_elevation = newz
-        ref_points = []
-        new_ref_points = []
-        x0 = 2002
-        y = 3007
-        yvec = range(4)
-        xvec = range(6)
-        for i in range(4):
-            y = y - 1
-            ynew = y - 3003.0
-            for j in range(6):
-                x = x0 + xvec[j]
-                xnew = x - 2002.0
-                if j <> 2 and (i<>1 or j<>4) and i<>3:
-                    ref_points.append([x,y])
-                    new_ref_points.append ([xnew,ynew])
-
-
-        #print points[:],points[:].shape
-        #print new_ref_points, len(new_ref_points)
-
-        assert num.allclose(elevation, ref_elevation)
-        assert num.allclose(points, new_ref_points)
-
-
-        #Cleanup
-        fid.close()
-
-
-        os.remove(root + '.pts')
-        os.remove(root + '.dem')
-        os.remove(root + '.asc')
-        os.remove(root + '.prj')
 
 
