Index: /trunk/anuga_core/source/anuga/file/sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/sww.py	(revision 7766)
+++ /trunk/anuga_core/source/anuga/file/sww.py	(revision 7767)
@@ -984,2 +984,328 @@
         log.critical('------------------------------------------------')
 
+
+
+
+##
+# @brief Get the extents of a NetCDF data file.
+# @param file_name The path to the SWW file.
+# @return A list of x, y, z and stage limits (min, max).
+def extent_sww(file_name):
+    """Read in an sww file.
+
+    Input:
+    file_name - the sww file
+
+    Output:
+    A list: [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)]
+    """
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    #Get NetCDF
+    fid = NetCDFFile(file_name, netcdf_mode_r)
+
+    # Get the variables
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    stage = fid.variables['stage'][:]
+
+    fid.close()
+
+    return [min(x), max(x), min(y), max(y), num.min(stage), num.max(stage)]
+
+
+##
+# @brief 
+# @param filename
+# @param boundary
+# @param t
+# @param fail_if_NaN
+# @param NaN_filler
+# @param verbose
+# @param very_verbose
+# @return 
+def load_sww_as_domain(filename, boundary=None, t=None,
+               fail_if_NaN=True, NaN_filler=0,
+               verbose=False, very_verbose=False):
+    """
+    Usage: domain = load_sww_as_domain('file.sww',
+                        t=time (default = last time in file))
+
+    Boundary is not recommended if domain.smooth is not selected, as it
+    uses unique coordinates, but not unique boundaries. This means that
+    the boundary file will not be compatable with the coordinates, and will
+    give a different final boundary, or crash.
+    """
+    
+    from Scientific.IO.NetCDF import NetCDFFile
+    from shallow_water import Domain
+
+    # initialise NaN.
+    NaN = 9.969209968386869e+036
+
+    if verbose: log.critical('Reading from %s' % filename)
+
+    fid = NetCDFFile(filename, netcdf_mode_r)    # Open existing file for read
+    time = fid.variables['time']       # Timesteps
+    if t is None:
+        t = time[-1]
+    time_interp = get_time_interp(time,t)
+
+    # Get the variables as numeric arrays
+    x = fid.variables['x'][:]                   # x-coordinates of vertices
+    y = fid.variables['y'][:]                   # y-coordinates of vertices
+    elevation = fid.variables['elevation']      # Elevation
+    stage = fid.variables['stage']              # Water level
+    xmomentum = fid.variables['xmomentum']      # Momentum in the x-direction
+    ymomentum = fid.variables['ymomentum']      # Momentum in the y-direction
+
+    starttime = fid.starttime[0]
+    volumes = fid.variables['volumes'][:]       # Connectivity
+    coordinates = num.transpose(num.asarray([x.tolist(), y.tolist()]))
+    # FIXME (Ole): Something like this might be better:
+    #                 concatenate((x, y), axis=1)
+    # or              concatenate((x[:,num.newaxis], x[:,num.newaxis]), axis=1)
+
+    conserved_quantities = []
+    interpolated_quantities = {}
+    other_quantities = []
+
+    # get geo_reference
+    try:                             # sww files don't have to have a geo_ref
+        geo_reference = Geo_reference(NetCDFObject=fid)
+    except: # AttributeError, e:
+        geo_reference = None
+
+    if verbose: log.critical('    getting quantities')
+
+    for quantity in fid.variables.keys():
+        dimensions = fid.variables[quantity].dimensions
+        if 'number_of_timesteps' in dimensions:
+            conserved_quantities.append(quantity)
+            interpolated_quantities[quantity] = \
+                  interpolated_quantity(fid.variables[quantity][:], time_interp)
+        else:
+            other_quantities.append(quantity)
+
+    other_quantities.remove('x')
+    other_quantities.remove('y')
+    #other_quantities.remove('z')
+    other_quantities.remove('volumes')
+    try:
+        other_quantities.remove('stage_range')
+        other_quantities.remove('xmomentum_range')
+        other_quantities.remove('ymomentum_range')
+        other_quantities.remove('elevation_range')
+    except:
+        pass
+
+    conserved_quantities.remove('time')
+
+    if verbose: log.critical('    building domain')
+
+    #    From domain.Domain:
+    #    domain = Domain(coordinates, volumes,\
+    #                    conserved_quantities = conserved_quantities,\
+    #                    other_quantities = other_quantities,zone=zone,\
+    #                    xllcorner=xllcorner, yllcorner=yllcorner)
+
+    # From shallow_water.Domain:
+    coordinates = coordinates.tolist()
+    volumes = volumes.tolist()
+    # FIXME:should this be in mesh? (peter row)
+    if fid.smoothing == 'Yes':
+        unique = False
+    else:
+        unique = True
+    if unique:
+        coordinates, volumes, boundary = weed(coordinates, volumes,boundary)
+
+      
+    
+    try:
+        domain = Domain(coordinates, volumes, boundary)
+    except AssertionError, e:
+        fid.close()
+        msg = 'Domain could not be created: %s. ' \
+              'Perhaps use "fail_if_NaN=False and NaN_filler = ..."' % e
+        raise DataDomainError, msg
+
+    if not boundary is None:
+        domain.boundary = boundary
+
+    domain.geo_reference = geo_reference
+
+    domain.starttime = float(starttime) + float(t)
+    domain.time = 0.0
+
+    for quantity in other_quantities:
+        try:
+            NaN = fid.variables[quantity].missing_value
+        except:
+            pass                       # quantity has no missing_value number
+        X = fid.variables[quantity][:]
+        if very_verbose:
+            log.critical('       %s' % str(quantity))
+            log.critical('        NaN = %s' % str(NaN))
+            log.critical('        max(X)')
+            log.critical('       %s' % str(max(X)))
+            log.critical('        max(X)==NaN')
+            log.critical('       %s' % str(max(X)==NaN))
+            log.critical('')
+        if max(X) == NaN or min(X) == NaN:
+            if fail_if_NaN:
+                msg = 'quantity "%s" contains no_data entry' % quantity
+                raise DataMissingValuesError, msg
+            else:
+                data = (X != NaN)
+                X = (X*data) + (data==0)*NaN_filler
+        if unique:
+            X = num.resize(X, (len(X)/3, 3))
+        domain.set_quantity(quantity, X)
+    #
+    for quantity in conserved_quantities:
+        try:
+            NaN = fid.variables[quantity].missing_value
+        except:
+            pass                       # quantity has no missing_value number
+        X = interpolated_quantities[quantity]
+        if very_verbose:
+            log.critical('       %s' % str(quantity))
+            log.critical('        NaN = %s' % str(NaN))
+            log.critical('        max(X)')
+            log.critical('       %s' % str(max(X)))
+            log.critical('        max(X)==NaN')
+            log.critical('       %s' % str(max(X)==NaN))
+            log.critical('')
+        if max(X) == NaN or min(X) == NaN:
+            if fail_if_NaN:
+                msg = 'quantity "%s" contains no_data entry' % quantity
+                raise DataMissingValuesError, msg
+            else:
+                data = (X != NaN)
+                X = (X*data) + (data==0)*NaN_filler
+        if unique:
+            X = num.resize(X, (X.shape[0]/3, 3))
+        domain.set_quantity(quantity, X)
+
+    fid.close()
+
+    return domain
+
+
+##
+# @brief 
+# @parm time 
+# @param t 
+# @return An (index, ration) tuple.
+def get_time_interp(time, t=None):
+    #Finds the ratio and index for time interpolation.
+    #It is borrowed from previous abstract_2d_finite_volumes code.
+    if t is None:
+        t=time[-1]
+        index = -1
+        ratio = 0.
+    else:
+        T = time
+        tau = t
+        index=0
+        msg = 'Time interval derived from file %s [%s:%s]' \
+              % ('FIXMEfilename', T[0], T[-1])
+        msg += ' does not match model time: %s' % tau
+        if tau < time[0]: raise DataTimeError, msg
+        if tau > time[-1]: raise DataTimeError, msg
+        while tau > time[index]: index += 1
+        while tau < time[index]: index -= 1
+        if tau == time[index]:
+            #Protect against case where tau == time[-1] (last time)
+            # - also works in general when tau == time[i]
+            ratio = 0
+        else:
+            #t is now between index and index+1
+            ratio = (tau - time[index])/(time[index+1] - time[index])
+
+    return (index, ratio)
+
+
+
+##
+# @brief Interpolate a quantity wrt time.
+# @param saved_quantity The quantity to interpolate.
+# @param time_interp (index, ratio)
+# @return A vector of interpolated values.
+def interpolated_quantity(saved_quantity, time_interp):
+    '''Given an index and ratio, interpolate quantity with respect to time.'''
+
+    index, ratio = time_interp
+
+    Q = saved_quantity
+
+    if ratio > 0:
+        q = (1-ratio)*Q[index] + ratio*Q[index+1]
+    else:
+        q = Q[index]
+
+    #Return vector of interpolated values
+    return q
+
+
+
+
+##
+# @brief 
+# @param coordinates 
+# @param volumes 
+# @param boundary 
+def weed(coordinates, volumes, boundary=None):
+    if isinstance(coordinates, num.ndarray):
+        coordinates = coordinates.tolist()
+    if isinstance(volumes, num.ndarray):
+        volumes = volumes.tolist()
+
+    unique = False
+    point_dict = {}
+    same_point = {}
+    for i in range(len(coordinates)):
+        point = tuple(coordinates[i])
+        if point_dict.has_key(point):
+            unique = True
+            same_point[i] = point
+            #to change all point i references to point j
+        else:
+            point_dict[point] = i
+            same_point[i] = point
+
+    coordinates = []
+    i = 0
+    for point in point_dict.keys():
+        point = tuple(point)
+        coordinates.append(list(point))
+        point_dict[point] = i
+        i += 1
+
+    for volume in volumes:
+        for i in range(len(volume)):
+            index = volume[i]
+            if index > -1:
+                volume[i] = point_dict[same_point[index]]
+
+    new_boundary = {}
+    if not boundary is None:
+        for segment in boundary.keys():
+            point0 = point_dict[same_point[segment[0]]]
+            point1 = point_dict[same_point[segment[1]]]
+            label = boundary[segment]
+            #FIXME should the bounday attributes be concaterated
+            #('exterior, pond') or replaced ('pond')(peter row)
+
+            if new_boundary.has_key((point0, point1)):
+                new_boundary[(point0,point1)] = new_boundary[(point0, point1)]
+
+            elif new_boundary.has_key((point1, point0)):
+                new_boundary[(point1,point0)] = new_boundary[(point1, point0)]
+            else: new_boundary[(point0, point1)] = label
+
+        boundary = new_boundary
+
+    return coordinates, volumes, boundary
Index: /trunk/anuga_core/source/anuga/file/test_sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7767)
+++ /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7767)
@@ -0,0 +1,165 @@
+import os
+import unittest
+import tempfile
+import numpy as num
+
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+from csv_file import load_csv_as_array, load_csv_as_dict
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
+from anuga.shallow_water.shallow_water_domain import Domain
+from sww import load_sww_as_domain
+
+# boundary functions
+from anuga.shallow_water.boundaries import Reflective_boundary, \
+            Field_boundary, Transmissive_momentum_set_stage_boundary, \
+            Transmissive_stage_zero_momentum_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import Transmissive_boundary, Dirichlet_boundary, \
+            Time_boundary, File_boundary, AWI_boundary
+
+
+class Test_sww(unittest.TestCase):
+    def setUp(self):
+        self.verbose = False
+        pass
+
+    def tearDown(self):
+        pass
+        
+    def test_sww2domain1(self):
+        ################################################
+        #Create a test domain, and evolve and save it.
+        ################################################
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+
+        yiel=0.01
+        points, vertices, boundary = rectangular(10,10)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.geo_reference = Geo_reference(56,11,11)
+        domain.smooth = False
+        domain.store = True
+        domain.set_name('bedslope')
+        domain.default_order=2
+        #Bed-slope and friction
+        domain.set_quantity('elevation', lambda x,y: -x/3)
+        domain.set_quantity('friction', 0.1)
+        # Boundary conditions
+        from math import sin, pi
+        Br = Reflective_boundary(domain)
+        Bt = Transmissive_boundary(domain)
+        Bd = Dirichlet_boundary([0.2,0.,0.])
+        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
+
+        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+
+        domain.quantities_to_be_stored['xmomentum'] = 2
+        domain.quantities_to_be_stored['ymomentum'] = 2
+        #Initial condition
+        h = 0.05
+        elevation = domain.quantities['elevation'].vertex_values
+        domain.set_quantity('stage', elevation + h)
+
+        domain.check_integrity()
+        #Evolution
+        #domain.tight_slope_limiters = 1
+        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
+            #domain.write_time()
+            pass
+
+
+
+        filename = domain.datadir + os.sep + domain.get_name() + '.sww'
+        domain2 = load_sww_as_domain(filename, None, fail_if_NaN=False,
+                                        verbose=self.verbose)
+        #points, vertices, boundary = rectangular(15,15)
+        #domain2.boundary = boundary
+        ###################
+        ##NOW TEST IT!!!
+        ###################
+
+        os.remove(filename)
+
+        bits = ['vertex_coordinates']
+        for quantity in domain.quantities_to_be_stored:
+            bits.append('get_quantity("%s").get_integral()' % quantity)
+            bits.append('get_quantity("%s").get_values()' % quantity)
+
+        for bit in bits:
+            #print 'testing that domain.'+bit+' has been restored'
+            #print bit
+            #print 'done'
+            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit))
+
+        ######################################
+        #Now evolve them both, just to be sure
+        ######################################x
+        domain.time = 0.
+        from time import sleep
+
+        final = .1
+        domain.set_quantity('friction', 0.1)
+        domain.store = False
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+
+
+        for t in domain.evolve(yieldstep = yiel, finaltime = final):
+            #domain.write_time()
+            pass
+
+        final = final - (domain2.starttime-domain.starttime)
+        #BUT since domain1 gets time hacked back to 0:
+        final = final + (domain2.starttime-domain.starttime)
+
+        domain2.smooth = False
+        domain2.store = False
+        domain2.default_order=2
+        domain2.set_quantity('friction', 0.1)
+        #Bed-slope and friction
+        # Boundary conditions
+        Bd2=Dirichlet_boundary([0.2,0.,0.])
+        domain2.boundary = domain.boundary
+        #print 'domain2.boundary'
+        #print domain2.boundary
+        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
+        #domain2.set_boundary({'exterior': Bd})
+
+        domain2.check_integrity()
+
+        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
+            #domain2.write_time()
+            pass
+
+        ###################
+        ##NOW TEST IT!!!
+        ##################
+
+        bits = ['vertex_coordinates']
+
+        for quantity in ['elevation','stage', 'ymomentum','xmomentum']:
+            bits.append('get_quantity("%s").get_integral()' %quantity)
+            bits.append('get_quantity("%s").get_values()' %quantity)
+
+        #print bits
+        for bit in bits:
+            #print bit
+            #print eval('domain.'+bit)
+            #print eval('domain2.'+bit)
+            
+            #print eval('domain.'+bit+'-domain2.'+bit)
+            msg = 'Values in the two domains are different for ' + bit
+            assert num.allclose(eval('domain.'+bit),eval('domain2.'+bit),
+                                rtol=1.e-5, atol=3.e-8), msg
+
+
+
+#################################################################################
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_sww, 'test')
+    runner = unittest.TextTestRunner(verbosity=1)
+    runner.run(suite)
