Index: inundation/pyvolution/data_manager.py
===================================================================
--- inundation/pyvolution/data_manager.py	(revision 2334)
+++ inundation/pyvolution/data_manager.py	(revision 2344)
@@ -2834,180 +2834,4 @@
 
 
-
-
-
-
-
-
-
-
-###############################################
-#OBSOLETE STUFF
-#Native checkpoint format.
-#Information needed to recreate a state is preserved
-#FIXME: Rethink and maybe use netcdf format
-def cpt_variable_writer(filename, t, v0, v1, v2):
-    """Store all conserved quantities to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'cpt', M, t)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(N):
-            fid.write('%.16e ' %v0[i,j])
-        for j in range(N):
-            fid.write('%.16e ' %v1[i,j])
-        for j in range(N):
-            fid.write('%.16e ' %v2[i,j])
-
-        fid.write('\n')
-    fid.close()
-
-
-def cpt_variable_reader(filename, t, v0, v1, v2):
-    """Store all conserved quantities to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'cpt', M, t)
-    #print 'Reading from %s' %FN
-
-    fid = open(FN)
-
-
-    for i in range(M):
-        values = fid.readline().split() #Get one line
-
-        for j in range(N):
-            v0[i,j] = float(values[j])
-            v1[i,j] = float(values[3+j])
-            v2[i,j] = float(values[6+j])
-
-    fid.close()
-
-def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
-    """Writes x,y,z,z,z coordinates of triangles constituting the bed
-    elevation.
-    FIXME: Not in use pt
-    """
-
-    M, N = v0.shape
-
-
-    print X0
-    import sys; sys.exit()
-    FN = create_filename(filename, 'cpt', M)
-    print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(2):
-            fid.write('%.16e ' %X0[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v0[i,j])       #z,z,z,
-
-        for j in range(2):
-            fid.write('%.16e ' %X1[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v1[i,j])
-
-        for j in range(2):
-            fid.write('%.16e ' %X2[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v2[i,j])
-
-        fid.write('\n')
-    fid.close()
-
-
-
-#Function for storing out to e.g. visualisation
-#FIXME: Do we want this?
-#FIXME: Not done yet for this version
-def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
-    """Writes x,y,z coordinates of triangles constituting the bed elevation.
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'dat', M)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(2):
-            fid.write('%f ' %X0[i,j])   #x, y
-        fid.write('%f ' %v0[i,0])       #z
-
-        for j in range(2):
-            fid.write('%f ' %X1[i,j])   #x, y
-        fid.write('%f ' %v1[i,0])       #z
-
-        for j in range(2):
-            fid.write('%f ' %X2[i,j])   #x, y
-        fid.write('%f ' %v2[i,0])       #z
-
-        fid.write('\n')
-    fid.close()
-
-
-
-def dat_variable_writer(filename, t, v0, v1, v2):
-    """Store water height to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'dat', M, t)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        fid.write('%.4f ' %v0[i,0])
-        fid.write('%.4f ' %v1[i,0])
-        fid.write('%.4f ' %v2[i,0])
-
-        fid.write('\n')
-    fid.close()
-
-
-def read_sww(filename):
-    """Read sww Net CDF file containing Shallow Water Wave simulation
-
-    The integer array volumes is of shape Nx3 where N is the number of
-    triangles in the mesh.
-
-    Each entry in volumes is an index into the x,y arrays (the location).
-
-    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
-    number_of_timesteps, number_of_points.
-
-    The momentum is not always stored.
-
-    """
-    from Scientific.IO.NetCDF import NetCDFFile
-    print 'Reading from ', filename
-    fid = NetCDFFile(filename, 'r')    #Open existing file for read
-#latitude, longitude
-    # Get the variables as Numeric arrays
-    x = fid.variables['x']             #x-coordinates of vertices
-    y = fid.variables['y']             #y-coordinates of vertices
-    z = fid.variables['elevation']     #Elevation
-    time = fid.variables['time']       #Timesteps
-    stage = fid.variables['stage']     #Water level
-    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
-    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
-
-    volumes = fid.variables['volumes'] #Connectivity
-
-    #FIXME (Ole): What is this?
-    #             Why isn't anything returned?
-    #             Where's the unit test?
-
 def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
                  verbose=False):
@@ -3138,223 +2962,4 @@
 
 
-
-def sww2asc_obsolete(basename_in, basename_out = None,
-            quantity = None,
-            timestep = None,
-            reduction = None,
-            cellsize = 10,
-            verbose = False,
-            origin = None):
-    """Read SWW file and convert to Digitial Elevation model format (.asc)
-
-    Example:
-
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    Also write accompanying file with same basename_in but extension .prj
-    used to fix the UTM zone, datum, false northings and eastings.
-
-    The prj format is assumed to be as
-
-    Projection    UTM
-    Zone          56
-    Datum         WGS84
-    Zunits        NO
-    Units         METERS
-    Spheroid      WGS84
-    Xshift        0.0000000000
-    Yshift        10000000.0000000000
-    Parameters
-
-
-    if quantity is given, out values from quantity otherwise default to
-    elevation
-
-    if timestep (an index) is given, output quantity at that timestep
-
-    if reduction is given use that to reduce quantity over all timesteps.
-
-    """
-    from Numeric import array, Float, concatenate, NewAxis, zeros,\
-         sometrue
-    from utilities.polygon import inside_polygon
-
-    #FIXME: Should be variable
-    datum = 'WGS84'
-    false_easting = 500000
-    false_northing = 10000000
-
-    if quantity is None:
-        quantity = 'elevation'
-
-    if reduction is None:
-        reduction = max
-
-    if basename_out is None:
-        basename_out = basename_in + '_%s' %quantity
-
-    swwfile = basename_in + '.sww'
-    ascfile = basename_out + '.asc'
-    prjfile = basename_out + '.prj'
-
-
-    if verbose: print 'Reading from %s' %swwfile
-    #Read sww file
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(swwfile)
-
-    #Get extent and reference
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    volumes = fid.variables['volumes'][:]
-
-    ymin = min(y); ymax = max(y)
-    xmin = min(x); xmax = max(x)
-
-    number_of_timesteps = fid.dimensions['number_of_timesteps']
-    number_of_points = fid.dimensions['number_of_points']
-    if origin is None:
-
-        #Get geo_reference
-        #sww files don't have to have a geo_ref
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference() #Default georef object
-
-        xllcorner = geo_reference.get_xllcorner()
-        yllcorner = geo_reference.get_yllcorner()
-        zone = geo_reference.get_zone()
-    else:
-        zone = origin[0]
-        xllcorner = origin[1]
-        yllcorner = origin[2]
-
-
-    #Get quantity and reduce if applicable
-    if verbose: print 'Reading quantity %s' %quantity
-
-    if quantity.lower() == 'depth':
-        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
-    else:
-        q = fid.variables[quantity][:]
-
-
-    if len(q.shape) == 2:
-        if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = zeros( number_of_points, Float )
-
-        for k in range(number_of_points):
-            q_reduced[k] = reduction( q[:,k] )
-
-        q = q_reduced
-
-    #Now q has dimension: number_of_points
-
-    #Create grid and update xll/yll corner and x,y
-    if verbose: print 'Creating grid'
-    ncols = int((xmax-xmin)/cellsize)+1
-    nrows = int((ymax-ymin)/cellsize)+1
-
-    newxllcorner = xmin+xllcorner
-    newyllcorner = ymin+yllcorner
-
-    x = x+xllcorner-newxllcorner
-    y = y+yllcorner-newyllcorner
-
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
-    assert len(vertex_points.shape) == 2
-
-
-    from Numeric import zeros, Float
-    grid_points = zeros ( (ncols*nrows, 2), Float )
-
-
-    for i in xrange(nrows):
-        yg = i*cellsize
-        for j in xrange(ncols):
-            xg = j*cellsize
-            k = i*ncols + j
-
-            grid_points[k,0] = xg
-            grid_points[k,1] = yg
-
-    #Interpolate
-    from least_squares import Interpolation
-
-
-    #FIXME: This should be done with precrop = True, otherwise it'll
-    #take forever. With expand_search  set to False, some grid points might
-    #miss out....
-
-    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
-                           precrop = False, expand_search = True,
-                           verbose = verbose)
-
-    #Interpolate using quantity values
-    if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q).flat
-
-    #Write
-    #Write prj file
-    if verbose: print 'Writing %s' %prjfile
-    prjid = open(prjfile, 'w')
-    prjid.write('Projection    %s\n' %'UTM')
-    prjid.write('Zone          %d\n' %zone)
-    prjid.write('Datum         %s\n' %datum)
-    prjid.write('Zunits        NO\n')
-    prjid.write('Units         METERS\n')
-    prjid.write('Spheroid      %s\n' %datum)
-    prjid.write('Xshift        %d\n' %false_easting)
-    prjid.write('Yshift        %d\n' %false_northing)
-    prjid.write('Parameters\n')
-    prjid.close()
-
-
-
-    if verbose: print 'Writing %s' %ascfile
-    NODATA_value = -9999
-
-    ascid = open(ascfile, 'w')
-
-    ascid.write('ncols         %d\n' %ncols)
-    ascid.write('nrows         %d\n' %nrows)
-    ascid.write('xllcorner     %d\n' %newxllcorner)
-    ascid.write('yllcorner     %d\n' %newyllcorner)
-    ascid.write('cellsize      %f\n' %cellsize)
-    ascid.write('NODATA_value  %d\n' %NODATA_value)
-
-
-    #Get bounding polygon from mesh
-    P = interp.mesh.get_boundary_polygon()
-    inside_indices = inside_polygon(grid_points, P)
-
-    for i in range(nrows):
-        if verbose and i%((nrows+10)/10)==0:
-            print 'Doing row %d of %d' %(i, nrows)
-
-        for j in range(ncols):
-            index = (nrows-i-1)*ncols+j
-
-            if sometrue(inside_indices == index):
-                ascid.write('%f ' %grid_values[index])
-            else:
-                ascid.write('%d ' %NODATA_value)
-
-        ascid.write('\n')
-
-    #Close
-    ascid.close()
-    fid.close()
-
-#********************
-#*** END OF OBSOLETE FUNCTIONS
-#***************
 
 
Index: inundation/pyvolution/data_manager_obsolete_stuff.py
===================================================================
--- inundation/pyvolution/data_manager_obsolete_stuff.py	(revision 2344)
+++ inundation/pyvolution/data_manager_obsolete_stuff.py	(revision 2344)
@@ -0,0 +1,395 @@
+
+
+
+
+
+
+
+###############################################
+#OBSOLETE STUFF
+#Native checkpoint format.
+#Information needed to recreate a state is preserved
+#FIXME: Rethink and maybe use netcdf format
+def cpt_variable_writer(filename, t, v0, v1, v2):
+    """Store all conserved quantities to file
+    """
+
+    M, N = v0.shape
+
+    FN = create_filename(filename, 'cpt', M, t)
+    #print 'Writing to %s' %FN
+
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(N):
+            fid.write('%.16e ' %v0[i,j])
+        for j in range(N):
+            fid.write('%.16e ' %v1[i,j])
+        for j in range(N):
+            fid.write('%.16e ' %v2[i,j])
+
+        fid.write('\n')
+    fid.close()
+
+
+def cpt_variable_reader(filename, t, v0, v1, v2):
+    """Store all conserved quantities to file
+    """
+
+    M, N = v0.shape
+
+    FN = create_filename(filename, 'cpt', M, t)
+    #print 'Reading from %s' %FN
+
+    fid = open(FN)
+
+
+    for i in range(M):
+        values = fid.readline().split() #Get one line
+
+        for j in range(N):
+            v0[i,j] = float(values[j])
+            v1[i,j] = float(values[3+j])
+            v2[i,j] = float(values[6+j])
+
+    fid.close()
+
+def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
+    """Writes x,y,z,z,z coordinates of triangles constituting the bed
+    elevation.
+    FIXME: Not in use pt
+    """
+
+    M, N = v0.shape
+
+
+    print X0
+    import sys; sys.exit()
+    FN = create_filename(filename, 'cpt', M)
+    print 'Writing to %s' %FN
+
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(2):
+            fid.write('%.16e ' %X0[i,j])   #x, y
+        for j in range(N):
+            fid.write('%.16e ' %v0[i,j])       #z,z,z,
+
+        for j in range(2):
+            fid.write('%.16e ' %X1[i,j])   #x, y
+        for j in range(N):
+            fid.write('%.16e ' %v1[i,j])
+
+        for j in range(2):
+            fid.write('%.16e ' %X2[i,j])   #x, y
+        for j in range(N):
+            fid.write('%.16e ' %v2[i,j])
+
+        fid.write('\n')
+    fid.close()
+
+
+
+#Function for storing out to e.g. visualisation
+#FIXME: Do we want this?
+#FIXME: Not done yet for this version
+def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
+    """Writes x,y,z coordinates of triangles constituting the bed elevation.
+    """
+
+    M, N = v0.shape
+
+    FN = create_filename(filename, 'dat', M)
+    #print 'Writing to %s' %FN
+
+    fid = open(FN, 'w')
+    for i in range(M):
+        for j in range(2):
+            fid.write('%f ' %X0[i,j])   #x, y
+        fid.write('%f ' %v0[i,0])       #z
+
+        for j in range(2):
+            fid.write('%f ' %X1[i,j])   #x, y
+        fid.write('%f ' %v1[i,0])       #z
+
+        for j in range(2):
+            fid.write('%f ' %X2[i,j])   #x, y
+        fid.write('%f ' %v2[i,0])       #z
+
+        fid.write('\n')
+    fid.close()
+
+
+
+def dat_variable_writer(filename, t, v0, v1, v2):
+    """Store water height to file
+    """
+
+    M, N = v0.shape
+
+    FN = create_filename(filename, 'dat', M, t)
+    #print 'Writing to %s' %FN
+
+    fid = open(FN, 'w')
+    for i in range(M):
+        fid.write('%.4f ' %v0[i,0])
+        fid.write('%.4f ' %v1[i,0])
+        fid.write('%.4f ' %v2[i,0])
+
+        fid.write('\n')
+    fid.close()
+
+
+def read_sww(filename):
+    """Read sww Net CDF file containing Shallow Water Wave simulation
+
+    The integer array volumes is of shape Nx3 where N is the number of
+    triangles in the mesh.
+
+    Each entry in volumes is an index into the x,y arrays (the location).
+
+    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
+    number_of_timesteps, number_of_points.
+
+    The momentum is not always stored.
+
+    """
+    from Scientific.IO.NetCDF import NetCDFFile
+    print 'Reading from ', filename
+    fid = NetCDFFile(filename, 'r')    #Open existing file for read
+#latitude, longitude
+    # Get the variables as Numeric arrays
+    x = fid.variables['x']             #x-coordinates of vertices
+    y = fid.variables['y']             #y-coordinates of vertices
+    z = fid.variables['elevation']     #Elevation
+    time = fid.variables['time']       #Timesteps
+    stage = fid.variables['stage']     #Water level
+    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
+    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
+
+    volumes = fid.variables['volumes'] #Connectivity
+
+    #FIXME (Ole): What is this?
+    #             Why isn't anything returned?
+    #             Where's the unit test?
+
+
+
+def sww2asc_obsolete(basename_in, basename_out = None,
+            quantity = None,
+            timestep = None,
+            reduction = None,
+            cellsize = 10,
+            verbose = False,
+            origin = None):
+    """Read SWW file and convert to Digitial Elevation model format (.asc)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Also write accompanying file with same basename_in but extension .prj
+    used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
+
+
+    if quantity is given, out values from quantity otherwise default to
+    elevation
+
+    if timestep (an index) is given, output quantity at that timestep
+
+    if reduction is given use that to reduce quantity over all timesteps.
+
+    """
+    from Numeric import array, Float, concatenate, NewAxis, zeros,\
+         sometrue
+    from utilities.polygon import inside_polygon
+
+    #FIXME: Should be variable
+    datum = 'WGS84'
+    false_easting = 500000
+    false_northing = 10000000
+
+    if quantity is None:
+        quantity = 'elevation'
+
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' %quantity
+
+    swwfile = basename_in + '.sww'
+    ascfile = basename_out + '.asc'
+    prjfile = basename_out + '.prj'
+
+
+    if verbose: print 'Reading from %s' %swwfile
+    #Read sww file
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    #Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    ymin = min(y); ymax = max(y)
+    xmin = min(x); xmax = max(x)
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+    if origin is None:
+
+        #Get geo_reference
+        #sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() #Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+
+    #Get quantity and reduce if applicable
+    if verbose: print 'Reading quantity %s' %quantity
+
+    if quantity.lower() == 'depth':
+        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
+    else:
+        q = fid.variables[quantity][:]
+
+
+    if len(q.shape) == 2:
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
+
+        for k in range(number_of_points):
+            q_reduced[k] = reduction( q[:,k] )
+
+        q = q_reduced
+
+    #Now q has dimension: number_of_points
+
+    #Create grid and update xll/yll corner and x,y
+    if verbose: print 'Creating grid'
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    newxllcorner = xmin+xllcorner
+    newyllcorner = ymin+yllcorner
+
+    x = x+xllcorner-newxllcorner
+    y = y+yllcorner-newyllcorner
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+
+    from Numeric import zeros, Float
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        yg = i*cellsize
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+
+
+    #FIXME: This should be done with precrop = True, otherwise it'll
+    #take forever. With expand_search  set to False, some grid points might
+    #miss out....
+
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False, expand_search = True,
+                           verbose = verbose)
+
+    #Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    grid_values = interp.interpolate(q).flat
+
+    #Write
+    #Write prj file
+    if verbose: print 'Writing %s' %prjfile
+    prjid = open(prjfile, 'w')
+    prjid.write('Projection    %s\n' %'UTM')
+    prjid.write('Zone          %d\n' %zone)
+    prjid.write('Datum         %s\n' %datum)
+    prjid.write('Zunits        NO\n')
+    prjid.write('Units         METERS\n')
+    prjid.write('Spheroid      %s\n' %datum)
+    prjid.write('Xshift        %d\n' %false_easting)
+    prjid.write('Yshift        %d\n' %false_northing)
+    prjid.write('Parameters\n')
+    prjid.close()
+
+
+
+    if verbose: print 'Writing %s' %ascfile
+    NODATA_value = -9999
+
+    ascid = open(ascfile, 'w')
+
+    ascid.write('ncols         %d\n' %ncols)
+    ascid.write('nrows         %d\n' %nrows)
+    ascid.write('xllcorner     %d\n' %newxllcorner)
+    ascid.write('yllcorner     %d\n' %newyllcorner)
+    ascid.write('cellsize      %f\n' %cellsize)
+    ascid.write('NODATA_value  %d\n' %NODATA_value)
+
+
+    #Get bounding polygon from mesh
+    P = interp.mesh.get_boundary_polygon()
+    inside_indices = inside_polygon(grid_points, P)
+
+    for i in range(nrows):
+        if verbose and i%((nrows+10)/10)==0:
+            print 'Doing row %d of %d' %(i, nrows)
+
+        for j in range(ncols):
+            index = (nrows-i-1)*ncols+j
+
+            if sometrue(inside_indices == index):
+                ascid.write('%f ' %grid_values[index])
+            else:
+                ascid.write('%d ' %NODATA_value)
+
+        ascid.write('\n')
+
+    #Close
+    ascid.close()
+    fid.close()
+
+#********************
+#*** END OF OBSOLETE FUNCTIONS
+#***************
