Changeset 7866

Timestamp:

Jun 22, 2010, 12:04:32 PM (13 years ago)

Author:

hudson

Message:

More swb tests passing. Cleaned up some pylint errors.

Location:

trunk/anuga_core/source/anuga

Files:

• __init__.py (modified) (1 diff)
• file/test_csv.py (modified) (1 diff)
• file/test_sww.py (modified) (1 diff)
• file_conversion/file_conversion.py (modified) (5 diffs)
• file_conversion/test_sww2dem.py (modified) (1 diff)
• file_conversion/test_urs2sts.py (modified) (1 diff)
• geometry/quad.py (modified) (3 diffs)
• geometry/test_geometry.py (modified) (5 diffs)
• geometry/test_polygon.py (modified) (13 diffs)
• pmesh/mesh_quadtree.py (modified) (2 diffs)
• shallow_water/test_data_manager.py (modified) (7 diffs)
• shallow_water_balanced/test_swb_basic.py (modified) (3 diffs)
• shallow_water_balanced/test_swb_boundary_condition.py (modified) (1 diff)
• shallow_water_balanced/test_swb_conservation.py (modified) (24 diffs)
• shallow_water_balanced/test_swb_cross_sections.py (modified) (4 diffs)
• shallow_water_balanced/test_swb_distribute.py (modified) (1 diff)
• shallow_water_balanced/test_swb_forcing_terms.py (modified) (36 diffs)
• shallow_water_balanced/test_swb_reporting.py (modified) (2 diffs)
• utilities/treenode.py (deleted)

trunk/anuga_core/source/anuga/__init__.py

@@ -93 +93 @@
 # Forcing
 #-----------------------------
-from anuga.shallow_water.forcing import Inflow
+from anuga.shallow_water.forcing import Inflow, Rainfall
 
 #-----------------------------

trunk/anuga_core/source/anuga/file/test_csv.py

@@ -220 +220 @@
         
 
+    
+    def test_csv2polygons(self):
+        """test loading of a csv polygon file.
+        """
+        
+        path = get_pathname_from_package('anuga.shallow_water')                
+        testfile = os.path.join(path, 'polygon_values_example.csv')                
+
+        polygons, values = load_csv_as_polygons(testfile, 
+                                        value_name='floors')
+
+        assert len(polygons) == 7, 'Must have seven polygons'
+        assert len(values) == 7, 'Must have seven values'
+            
+        # Known floor values
+        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
+        
+        # Known polygon values
+        known_polys = {'1': [[422681.61,871117.55],
+                             [422691.02,871117.60],
+                             [422690.87,871084.23],
+                             [422649.36,871081.85],
+                             [422649.36,871080.39],
+                             [422631.86,871079.50],
+                             [422631.72,871086.75],
+                             [422636.75,871087.20],
+                             [422636.75,871091.50],
+                             [422649.66,871092.09],
+                             [422649.83,871084.91],
+                             [422652.94,871084.90],
+                             [422652.84,871092.39],
+                             [422681.83,871093.73],
+                             [422681.61,871117.55]],
+                       '2': [[422664.22,870785.46],
+                             [422672.48,870780.14],
+                             [422668.17,870772.62],
+                             [422660.35,870777.17],
+                             [422664.22,870785.46]],
+                       '3': [[422661.30,871215.06],
+                             [422667.50,871215.70],
+                             [422668.30,871204.86],
+                             [422662.21,871204.33],
+                             [422661.30,871215.06]],
+                       '4': [[422473.44,871191.22],
+                             [422478.33,871192.26],
+                             [422479.52,871186.03],
+                             [422474.78,871185.14],
+                             [422473.44,871191.22]],
+                       '5': [[422369.69,871049.29],
+                             [422378.63,871053.58],
+                             [422383.91,871044.51],
+                             [422374.97,871040.32],
+                             [422369.69,871049.29]],
+                       '8': [[422730.56,871203.13],
+                             [422734.10,871204.90],
+                             [422735.26,871202.18],
+                             [422731.87,871200.58],
+                             [422730.56,871203.13]],
+                       '9': [[422659.85,871213.80],
+                             [422660.91,871210.97],
+                             [422655.42,871208.85],
+                             [422654.36,871211.68],
+                             [422659.85,871213.80]]
+                       }        
+        
+
+        
+                
+        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
+            assert id in polygons.keys()
+            assert id in values.keys()            
+
+            assert int(values[id]) == int(floors[id])
+            assert len(polygons[id]) == len(known_polys[id])
+            assert num.allclose(polygons[id], known_polys[id])
+
+
+    def test_csv2polygons_with_clipping(self):
+        """test_csv2polygons with optional clipping
+        """
+        #FIXME(Ole): Not Done!!
+        
+        path = get_pathname_from_package('anuga.shallow_water')                
+        testfile = os.path.join(path, 'polygon_values_example.csv')                
+
+        polygons, values = load_csv_as_polygons(testfile, 
+                                        value_name='floors',
+                                        clipping_polygons=None)
+
+        assert len(polygons) == 7, 'Must have seven polygons'
+        assert len(values) == 7, 'Must have seven values'
+            
+        # Known floor values
+        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
+        
+        # Known polygon values
+        known_polys = {'1': [[422681.61,871117.55],
+                             [422691.02,871117.60],
+                             [422690.87,871084.23],
+                             [422649.36,871081.85],
+                             [422649.36,871080.39],
+                             [422631.86,871079.50],
+                             [422631.72,871086.75],
+                             [422636.75,871087.20],
+                             [422636.75,871091.50],
+                             [422649.66,871092.09],
+                             [422649.83,871084.91],
+                             [422652.94,871084.90],
+                             [422652.84,871092.39],
+                             [422681.83,871093.73],
+                             [422681.61,871117.55]],
+                       '2': [[422664.22,870785.46],
+                             [422672.48,870780.14],
+                             [422668.17,870772.62],
+                             [422660.35,870777.17],
+                             [422664.22,870785.46]],
+                       '3': [[422661.30,871215.06],
+                             [422667.50,871215.70],
+                             [422668.30,871204.86],
+                             [422662.21,871204.33],
+                             [422661.30,871215.06]],
+                       '4': [[422473.44,871191.22],
+                             [422478.33,871192.26],
+                             [422479.52,871186.03],
+                             [422474.78,871185.14],
+                             [422473.44,871191.22]],
+                       '5': [[422369.69,871049.29],
+                             [422378.63,871053.58],
+                             [422383.91,871044.51],
+                             [422374.97,871040.32],
+                             [422369.69,871049.29]],
+                       '8': [[422730.56,871203.13],
+                             [422734.10,871204.90],
+                             [422735.26,871202.18],
+                             [422731.87,871200.58],
+                             [422730.56,871203.13]],
+                       '9': [[422659.85,871213.80],
+                             [422660.91,871210.97],
+                             [422655.42,871208.85],
+                             [422654.36,871211.68],
+                             [422659.85,871213.80]]
+                       }        
+        
+
+        
+                
+        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
+            assert id in polygons.keys()
+            assert id in values.keys()            
+
+            assert int(values[id]) == int(floors[id])
+            assert len(polygons[id]) == len(known_polys[id])
+            assert num.allclose(polygons[id], known_polys[id])
+
+
+
+
+    
+    def test_csv2building_polygons(self):
+        """test_csv2building_polygons
+        """
+        
+        path = get_pathname_from_package('anuga.shallow_water')                
+        testfile = os.path.join(path, 'polygon_values_example.csv')                
+
+        polygons, values = load_csv_as_building_polygons(testfile, 
+                                                 floor_height=3)
+
+        assert len(polygons) == 7, 'Must have seven polygons'
+        assert len(values) == 7, 'Must have seven values'
+            
+        # Known floor values
+        floors = {'1': 6, '2': 0, '3': 3, '4': 6, '5': 0, '8': 3, '9': 3}
+        
+                
+        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
+            assert id in polygons.keys()
+            assert id in values.keys()            
+
+            assert float(values[id]) == float(floors[id])
+
 
 

trunk/anuga_core/source/anuga/file/test_sww.py

@@ -251 +251 @@
                 assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
 
+
+    def test_triangulation(self):
+        # 
+        #  
+        
+        filename = tempfile.mktemp("_data_manager.sww")
+        outfile = NetCDFFile(filename, netcdf_mode_w)
+        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
+        volumes = (0,1,2)
+        elevation = [0,1,2]
+        new_origin = None
+        new_origin = Geo_reference(56, 0, 0)
+        times = [0, 10]
+        number_of_volumes = len(volumes)
+        number_of_points = len(points_utm)
+        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
+        sww.store_header(outfile, times, number_of_volumes,
+                         number_of_points, description='fully sick testing',
+                         verbose=self.verbose,sww_precision=netcdf_float)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                verbose=self.verbose)       
+        outfile.close()
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        fid.close()
+
+        assert num.allclose(num.array(map(None, x,y)), points_utm)
+        os.remove(filename)
+
+        
+    def test_triangulationII(self):
+        # 
+        #  
+        
+        filename = tempfile.mktemp("_data_manager.sww")
+        outfile = NetCDFFile(filename, netcdf_mode_w)
+        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
+        volumes = (0,1,2)
+        elevation = [0,1,2]
+        new_origin = None
+        #new_origin = Geo_reference(56, 0, 0)
+        times = [0, 10]
+        number_of_volumes = len(volumes)
+        number_of_points = len(points_utm)
+        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
+        sww.store_header(outfile, times, number_of_volumes,
+                         number_of_points, description='fully sick testing',
+                         verbose=self.verbose,sww_precision=netcdf_float)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                new_origin=new_origin,
+                                verbose=self.verbose)
+        sww.store_static_quantities(outfile, elevation=elevation)                                
+                                
+        outfile.close()
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        results_georef = Geo_reference()
+        results_georef.read_NetCDF(fid)
+        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
+        fid.close()
+
+        assert num.allclose(num.array(map(None, x,y)), points_utm)
+        os.remove(filename)
+
+        
+    def test_triangulation_new_origin(self):
+        # 
+        #  
+        
+        filename = tempfile.mktemp("_data_manager.sww")
+        outfile = NetCDFFile(filename, netcdf_mode_w)
+        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
+        volumes = (0,1,2)
+        elevation = [0,1,2]
+        new_origin = None
+        new_origin = Geo_reference(56, 1, 554354)
+        points_utm = new_origin.change_points_geo_ref(points_utm)
+        times = [0, 10]
+        number_of_volumes = len(volumes)
+        number_of_points = len(points_utm)
+        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
+        sww.store_header(outfile, times, number_of_volumes,
+                         number_of_points, description='fully sick testing',
+                         verbose=self.verbose,sww_precision=netcdf_float)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                verbose=self.verbose)
+        outfile.close()
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        results_georef = Geo_reference()
+        results_georef.read_NetCDF(fid)
+        assert results_georef == new_origin
+        fid.close()
+
+        absolute = Geo_reference(56, 0,0)
+        assert num.allclose(num.array(
+            absolute.change_points_geo_ref(map(None, x,y),
+                                           new_origin)),points_utm)
+        
+        os.remove(filename)
+        
+    def test_triangulation_points_georeference(self):
+        # 
+        #  
+        
+        filename = tempfile.mktemp("_data_manager.sww")
+        outfile = NetCDFFile(filename, netcdf_mode_w)
+        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
+        volumes = (0,1,2)
+        elevation = [0,1,2]
+        new_origin = None
+        points_georeference = Geo_reference(56, 1, 554354)
+        points_utm = points_georeference.change_points_geo_ref(points_utm)
+        times = [0, 10]
+        number_of_volumes = len(volumes)
+        number_of_points = len(points_utm)
+        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
+        sww.store_header(outfile, times, number_of_volumes,
+                         number_of_points, description='fully sick testing',
+                         verbose=self.verbose,sww_precision=netcdf_float)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                points_georeference=points_georeference,
+                                verbose=self.verbose)       
+        outfile.close()
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        results_georef = Geo_reference()
+        results_georef.read_NetCDF(fid)
+        assert results_georef == points_georeference
+        fid.close()
+
+        assert num.allclose(num.array(map(None, x,y)), points_utm)
+        os.remove(filename)
+        
+    def test_triangulation_2_geo_refs(self):
+        # 
+        #  
+        
+        filename = tempfile.mktemp("_data_manager.sww")
+        outfile = NetCDFFile(filename, netcdf_mode_w)
+        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
+        volumes = (0,1,2)
+        elevation = [0,1,2]
+        new_origin = Geo_reference(56, 1, 1)
+        points_georeference = Geo_reference(56, 0, 0)
+        points_utm = points_georeference.change_points_geo_ref(points_utm)
+        times = [0, 10]
+        number_of_volumes = len(volumes)
+        number_of_points = len(points_utm)
+        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
+        sww.store_header(outfile, times, number_of_volumes,
+                         number_of_points, description='fully sick testing',
+                         verbose=self.verbose,sww_precision=netcdf_float)
+        sww.store_triangulation(outfile, points_utm, volumes,
+                                elevation,  new_origin=new_origin,
+                                points_georeference=points_georeference,
+                                verbose=self.verbose)       
+        outfile.close()
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        results_georef = Geo_reference()
+        results_georef.read_NetCDF(fid)
+        assert results_georef == new_origin
+        fid.close()
+
+
+        absolute = Geo_reference(56, 0,0)
+        assert num.allclose(num.array(
+            absolute.change_points_geo_ref(map(None, x,y),
+                                           new_origin)),points_utm)
+        os.remove(filename)
+
 #################################################################################
 

trunk/anuga_core/source/anuga/file_conversion/file_conversion.py

@@ -91 +91 @@
 # @param quantity_names 
 # @param time_as_seconds 
-def timefile2netcdf(file_text, quantity_names=None, time_as_seconds=False):
+def timefile2netcdf(file_text, file_out = None, quantity_names=None, \
+                                time_as_seconds=False):
     """Template for converting typical text files with time series to
     NetCDF tms file.
@@ -111 +112 @@
     will provide a time dependent function f(t) with three attributes
 
-    filename is assumed to be the rootname with extenisons .txt and .sww
+    filename is assumed to be the rootname with extensions .txt/.tms and .sww
     """
 
@@ -121 +122 @@
         raise IOError('Input file %s should be of type .txt.' % file_text)
 
-    filename = file_text[:-4]
+    if file_out == None:
+        file_out = file_text[:-4] + '.tms'
+
     fid = open(file_text)
     line = fid.readline()
@@ -181 +184 @@
             Q[i, j] = float(value)
 
-    msg = 'File %s must list time as a monotonuosly ' % filename
+    msg = 'File %s must list time as a monotonuosly ' % file_text
     msg += 'increasing sequence'
     assert num.alltrue(T[1:] - T[:-1] > 0), msg
@@ -188 +191 @@
     from Scientific.IO.NetCDF import NetCDFFile
 
-    fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
+    fid = NetCDFFile(file_out, netcdf_mode_w)
 
     fid.institution = 'Geoscience Australia'

trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py

@@ -1520 +1520 @@
         os.remove(ascfile)
         os.remove(swwfile2)
+
+
+    def test_export_grid(self):
+        """
+        test_export_grid(self):
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        self.domain.set_name('teg')
+
+        prjfile = self.domain.get_name() + '_elevation.prj'
+        ascfile = self.domain.get_name() + '_elevation.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        fid.close()
+
+        #Export to ascii/prj files
+        sww2dem_batch(self.domain.get_name(),
+                quantities = 'elevation',
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_export_gridII(self):
+        """
+        test_export_gridII(self):
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        self.domain.set_name('tegII')
+
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]        
+
+        #print 'stage', stage
+        #print 'xmom', xmomentum
+        #print 'ymom', ymomentum        
+
+        fid.close()
+
+        #Export to ascii/prj files
+        if True:
+            sww2dem_batch(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['depth'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+
+            export_grid(self.domain.get_name(),
+                quantities = ['elevation'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        prjfile = self.domain.get_name() + '_elevation.prj'
+        ascfile = self.domain.get_name() + '_elevation.asc'
+        
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #print "ascfile", ascfile
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        
+        #Check asc file
+        ascfile = self.domain.get_name() + '_depth.asc'
+        prjfile = self.domain.get_name() + '_depth.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])                
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+    def test_export_gridIII(self):
+        """
+        test_export_gridIII
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        
+        self.domain.set_name('tegIII')
+
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+        
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep() #'stage')
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep() #'stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        fid.close()
+
+        #Export to ascii/prj files
+        extra_name_out = 'yeah'
+        if True:
+            sww2dem_batch(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        extra_name_out = extra_name_out,
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['depth'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['elevation'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
+        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
+        
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #print "ascfile", ascfile
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        
+        #Check asc file
+        ascfile = self.domain.get_name() + '_depth_yeah.asc'
+        prjfile = self.domain.get_name() + '_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_export_grid_bad(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        try:
+            sww2dem_batch('a_small_round-egg',
+                        quantities = ['elevation', 'depth'],
+                        cellsize = 99,
+                        verbose = self.verbose,
+                        format = 'asc')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
         
 

trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py

@@ -1958 +1958 @@
         os.remove(meshname)
         
-
+        
+        
+    def test_file_boundary_sts_time_limit(self):
+        """test_file_boundary_stsIV_sinewave_ordering(self):
+        Read correct points from ordering file and apply sts to boundary
+        This one uses a sine wave and compares to time boundary
+        
+        This one tests that times used can be limited by upper limit
+        """
+        
+        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide = 0.35
+        time_step_count = 50
+        time_step = 0.1
+        times_ref = num.arange(0, time_step_count*time_step, time_step)
+        
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        
+        gauge_depth=20*num.ones(n,num.float)
+        
+        ha1=num.ones((n,time_step_count),num.float)
+        ua1=3.*num.ones((n,time_step_count),num.float)
+        va1=2.*num.ones((n,time_step_count),num.float)
+        for i in range(n):
+            ha1[i]=num.sin(times_ref)
+        
+        
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha1,
+                                           ua=ua1,
+                                           va=va1)
+
+        # Write order file
+        file_handle, order_base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(order_base_name)
+        d=","
+        order_file=order_base_name+'order.txt'
+        fid=open(order_file,'w')
+        
+        # Write Header
+        header='index, longitude, latitude\n'
+        fid.write(header)
+        indices=[3,0,1]
+        for i in indices:
+            line=str(i)+d+str(lat_long_points[i][1])+d+\
+                str(lat_long_points[i][0])+"\n"
+            fid.write(line)
+        fid.close()
+
+        sts_file=base_name
+        urs2sts(base_name, basename_out=sts_file,
+                ordering_filename=order_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+        
+        
+        
+        # Now read the sts file and check that values have been stored correctly.
+        fid = NetCDFFile(sts_file + '.sts')
+
+        # Check the time vector
+        times = fid.variables['time'][:]
+        starttime = fid.starttime[0]
+        #print times
+        #print starttime
+
+        # Check sts quantities
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        
+
+        # Create beginnings of boundary polygon based on sts_boundary
+        boundary_polygon = create_sts_boundary(base_name)
+        
+        os.remove(order_file)
+
+        # Append the remaining part of the boundary polygon to be defined by
+        # the user
+        bounding_polygon_utm=[]
+        for point in bounding_polygon:
+            zone,easting,northing=redfearn(point[0],point[1])
+            bounding_polygon_utm.append([easting,northing])
+
+        boundary_polygon.append(bounding_polygon_utm[3])
+        boundary_polygon.append(bounding_polygon_utm[4])
+
+        #print 'boundary_polygon', boundary_polygon
+        
+
+        assert num.allclose(bounding_polygon_utm,boundary_polygon)
+
+
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        
+        # have to change boundary tags from last example because now bounding
+        # polygon starts in different place.
+        create_mesh_from_regions(boundary_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=boundary_polygon)
+        time_vec = Bf.F.get_time()
+        assert num.allclose(Bf.F.starttime, starttime)
+        assert num.allclose(time_vec, times_ref)                                   
+        
+        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=time_limit+starttime,
+                               boundary_polygon=boundary_polygon)
+        
+            time_vec = Bf.F.get_time()
+            assert num.allclose(Bf.F.starttime, starttime)            
+            assert num.alltrue(time_vec < time_limit)
+            
+            
+        try:    
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=-1+starttime,
+                               boundary_polygon=boundary_polygon)            
+            time_vec = Bf.F.get_time()    
+            print time_vec    
+        except AssertionError:
+            pass
+        else:
+            raise Exception, 'Should have raised Exception here'
 
 #-------------------------------------------------------------

trunk/anuga_core/source/anuga/geometry/quad.py

@@ -13 +13 @@
 """
 
-from anuga.utilities.treenode import TreeNode
 import anuga.utilities.log as log
 
             
-class Cell(TreeNode):
+class Cell():
     """class Cell
 
@@ -29 +28 @@
             parent is the node above this one, or None if it is root.
         """
-  
-        # Initialise base classes
-        TreeNode.__init__(self, name)
     
         self.extents = extents
@@ -47 +43 @@
             ret_str += ', children: %d' % (len(self.children))
         return ret_str
-
-   
-
-    def clear(self):
-        """ Remove all leaves from this node.
-        """
-        self.Prune()   # TreeNode method
-
-
-    def clear_leaf_node(self):
-        """ Clears storage in leaf node.
-            Called from Treenode.
-            Must exist.    
-        """
-        self.leaves = []
-    
-    
-    def clear_internal_node(self):
-        """Called from Treenode.    
-    Must exist.
-    """
-        self.leaves = []
 
 

trunk/anuga_core/source/anuga/geometry/test_geometry.py

@@ +1 @@
+""" Test for the geometry classes.
+
+    Pylint quality rating as of June 2010: 8.51/10.
+"""
+
 import unittest
-import numpy as num
 
 from aabb import AABB
 from quad import Cell
 
-import types, sys
+import types
 
 #-------------------------------------------------------------
 
 class Test_Geometry(unittest.TestCase):
-
+    """ Test geometry classes. """
     def setUp(self):
+        """ Generic set up for geometry tests. """
         pass
 
     def tearDown(self):
+        """ Generic shut down for geometry tests. """
         pass
 
-    def test_AABB_contains(self):
+    def test_aabb_contains(self):
+        """ Test if point is correctly classified as falling inside or
+            outside of bounding box. """
         box = AABB(1, 21, 1, 11)
         assert box.contains([10, 5])
@@ -28 +36 @@
         assert not box.contains([50, 6])        
         
-    def test_AABB_split_vert(self):
+    def test_aabb_split_vert(self):
+        """ Test that a bounding box can be split correctly along an axis.
+        """
         parent = AABB(1, 21, 1, 11)
         
@@ -43 +53 @@
         self.assertEqual(child2.ymax, 11)    
 
-    def test_AABB_split_horiz(self):
+    def test_aabb_split_horiz(self):
+        """ Test that a bounding box will be split along the horizontal axis
+        correctly. """
         parent = AABB(1, 11, 1, 41)
         
@@ -59 +71 @@
         
     def test_add_data(self):
-        cell = Cell(AABB(0,10, 0,5), None)
-        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
+        """ Test add and retrieve arbitrary data from tree structure. """
+        cell = Cell(AABB(0, 10, 0, 5), None)
+        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), 222),  \
                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
 
         result = cell.retrieve()
-        assert type(result) in [types.ListType,types.TupleType],\
+        assert type(result) in [types.ListType,types.TupleType], \
                             'should be a list'
 
-        self.assertEqual(len(result),4)
+        self.assertEqual(len(result), 4)
         
     def test_search(self):
+        """ Test search tree for an intersection. """
         test_tag = 222
-        cell = Cell(AABB(0,10, 0,5), None)
-        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), test_tag),  \
+        cell = Cell(AABB(0, 10, 0,5), None)
+        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), test_tag),  \
                      (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
 
         result = cell.search([8.5, 1.5])
-        assert type(result) in [types.ListType,types.TupleType],\
+        assert type(result) in [types.ListType, types.TupleType], \
                             'should be a list'
         assert(len(result) == 1)
-        data, node = result[0]
+        data, _ = result[0]
         self.assertEqual(data, test_tag, 'only 1 point should intersect')
 
     def test_get_siblings(self):
         """ Make sure children know their parent. """
-        cell = Cell(AABB(0,10, 0,5), None)
-        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222)])
+        cell = Cell(AABB(0, 10, 0, 5), None)
+        cell.insert([(AABB(1, 3, 1, 3), 111), (AABB(8, 9, 1, 2), 222)])
 
         assert len(cell.children) == 2
@@ -92 +106 @@
         assert cell.children[1].parent == cell
 
-    def test_clear_1(self):
-        cell = Cell(AABB(0,10, 0,5), None)    
-        cell.insert([(AABB(1,3, 1, 3), 111), (AABB(8,9, 1, 2), 222),  \
-                     (AABB(7, 8, 3, 4), 333), (AABB(1, 10, 0, 1), 444)])
-                     
-        assert len(cell.retrieve()) == 4
-        cell.clear()
-
-        assert len(cell.retrieve()) == 0
 
 ################################################################################

trunk/anuga_core/source/anuga/geometry/test_polygon.py

@@ -1 +1 @@
 #!/usr/bin/env python
+
+""" Test suite to test polygon functionality. """
 
 import unittest
 import numpy as num
-from math import sqrt, pi
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.utilities.system_tools import get_pathname_from_package
 
-from polygon import *
-from polygon import _poly_xy
+from polygon import _poly_xy, separate_points_by_polygon, \
+                    populate_polygon, polygon_area, is_inside_polygon, \
+                    read_polygon, point_on_line, point_in_polygon, \
+                    plot_polygons, outside_polygon, is_outside_polygon, \
+                    intersection, is_complex, is_inside_triangle, \
+                    interpolate_polyline, inside_polygon, in_and_outside_polygon
+                    
 from polygon_function import Polygon_function
 from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -401 +407 @@
 
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], \
+                        [0.5, -0.5]]
         res, count = separate_points_by_polygon( points, polygon )
 
@@ -1401 +1408 @@
 
     def test_no_intersection(self):
+        """ Test 2 non-touching lines don't intersect. """
         line0 = [[-1,1], [1,1]]
         line1 = [[0,-1], [0,0]]
@@ -1469 +1477 @@
 
         # Overlap
-        line0 = [[0,5], [7,19]]
-        line1 = [[1,7], [10,25]]
+        line0 = [[0, 5], [7, 19]]
+        line1 = [[1, 7], [10, 25]]
         status, value = intersection(line0, line1)
         assert status == 2
@@ -1512 +1520 @@
 
 
-    def NOtest_inside_polygon_main(self):
-        # FIXME (Ole): Why is this disabled?
-        #print "inside",inside
-        #print "outside",outside
-
-        assert not inside_polygon((0.5, 1.5), polygon)
-        assert not inside_polygon((0.5, -0.5), polygon)
-        assert not inside_polygon((-0.5, 0.5), polygon)
-        assert not inside_polygon((1.5, 0.5), polygon)
-
-        # Try point on borders
-        assert inside_polygon((1., 0.5), polygon, closed=True)
-        assert inside_polygon((0.5, 1), polygon, closed=True)
-        assert inside_polygon((0., 0.5), polygon, closed=True)
-        assert inside_polygon((0.5, 0.), polygon, closed=True)
-
-        assert not inside_polygon((0.5, 1), polygon, closed=False)
-        assert not inside_polygon((0., 0.5), polygon, closed=False)
-        assert not inside_polygon((0.5, 0.), polygon, closed=False)
-        assert not inside_polygon((1., 0.5), polygon, closed=False)
-
+    def test_inside_polygon_main(self):
         # From real example (that failed)
         polygon = [[20,20], [40,20], [40,40], [20,40]]
@@ -1546 +1534 @@
 
     def test_polygon_area(self):
+        """ Test getting the area of a polygon. """
         # Simplest case: Polygon is the unit square
         polygon = [[0,0], [1,0], [1,1], [0,1]]
@@ -1611 +1600 @@
         # Another case
         polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
-        v = plot_polygons([polygon2,polygon3], figname='test2')
+        plot_polygons([polygon2, polygon3], figname='test2')
 
         for file in ['test1.png', 'test2.png']:
@@ -1619 +1608 @@
 
     def test_inside_polygon_geospatial(self):
+        """ Test geospatial coords inside poly. """
         #Simplest case: Polygon is the unit square
-        polygon_absolute = [[0,0], [1,0], [1,1], [0,1]]
+        polygon_absolute = [[0, 0], [1, 0], [1, 1], [0, 1]]
         poly_geo_ref = Geo_reference(57, 100, 100)
         polygon = poly_geo_ref.change_points_geo_ref(polygon_absolute)
@@ -1639 +1629 @@
         assert is_inside_polygon(points_absolute, polygon_absolute)
 
-    def NOtest_decimate_polygon(self):
+    def test_decimate_polygon(self):
+        from polygon import decimate_polygon
         polygon = [[0,0], [10,10], [15,5], [20, 10],
                    [25,0], [30,10], [40,-10], [35, -5]]
@@ -1645 +1636 @@
         dpoly = decimate_polygon(polygon, factor=2)
 
-        print dpoly
-
         assert len(dpoly)*2==len(polygon)
 
-        minx = maxx = polygon[0][0]
-        miny = maxy = polygon[0][1]
-        for point in polygon[1:]:
-            x, y = point
-
-            if x < minx: minx = x
-            if x > maxx: maxx = x
-            if y < miny: miny = y
-            if y > maxy: maxy = y
-
-        assert [minx, miny] in polygon
-        print minx, maxy
-        assert [minx, maxy] in polygon
-        assert [maxx, miny] in polygon
-        assert [maxx, maxy] in polygon
 
     def test_interpolate_polyline(self):
@@ -1756 +1730 @@
 
     def test_is_inside_triangle_more(self):
-
+        """ Test if points inside triangles are detected correctly. """
         res = is_inside_triangle([0.5, 0.5], [[ 0.5,  0. ],
                                               [ 0.5,  0.5],
@@ -1802 +1776 @@
         
     def test_is_polygon_complex(self):
-        """Test a concave and a complex poly with is_complex, to make sure it can detect
-           self-intersection.
+        """ Test a concave and a complex poly with is_complex, to make
+            sure it can detect self-intersection.
         """
         concave_poly = [[0, 0], [10, 0], [5, 5], [10, 10], [0, 10]]
-        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], [0, 10]]
+        complex_poly = [[0, 0], [10, 0], [5, 5], [4, 15], [5, 7], [10, 10], \
+                            [0, 10]]
 
         assert not is_complex(concave_poly)
@@ -1812 +1787 @@
 
     def test_is_polygon_complex2(self):
-        """Test a concave and a complex poly with is_complex, to make sure it can detect
-           self-intersection. This test uses more complicated polygons.
+        """ Test a concave and a complex poly with is_complex, to make sure it
+            can detect self-intersection. This test uses more complicated
+            polygons.
         """    
-        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7,6], [10, 10], [1,5], [0, 10]]
-        complex_poly = [[0, 0], [12,12], [10, 0], [5, 5], [3,18], [4, 15], [5, 7], [10, 10], [0, 10], [16, 12]]       
+        concave_poly = [[0, 0], [10, 0], [11,0], [5, 5], [7, 6], [10, 10], \
+                        [1, 5], [0, 10]]
+        complex_poly = [[0, 0], [12, 12], [10, 0], [5, 5], [3,18], [4, 15], \
+                        [5, 7], [10, 10], [0, 10], [16, 12]]       
 
         assert not is_complex(concave_poly)

trunk/anuga_core/source/anuga/pmesh/mesh_quadtree.py

@@ -18 +18 @@
     from polygon_ext import _is_inside_triangle        
 else:
-    msg = 'C implementations could not be accessed by %s.\n ' % __file__
-    msg += 'Make sure compile_all.py has been run as described in '
-    msg += 'the ANUGA installation guide.'
-    raise Exception, msg
+    MESSAGE = 'C implementations could not be accessed by %s.\n ' % __file__
+    MESSAGE += 'Make sure compile_all.py has been run as described in '
+    MESSAGE += 'the ANUGA installation guide.'
+    raise Exception(MESSAGE)
 
 import numpy as num
@@ -50 +50 @@
         extents.grow(1.001) # To avoid round off error
         Cell.__init__(self, extents, None)  # root has no parent
-        
+
+        self.last_triangle = None        
         N = len(mesh)
         self.mesh = mesh

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

@@ -503 +503 @@
         sww.store_timestep()
 
-
         # Check contents
         # Get NetCDF
@@ -530 +529 @@
         assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
                                    Q0[5] + Q2[6] + Q1[7])/6)
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(sww.filename)
-
-    def no_test_sww_variable3(self):
-        """Test that sww information can be written correctly
-        multiple timesteps using a different reduction operator (min)
-        """
-
-        # Different reduction in sww files has been made obsolete.
-        
-        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 = min
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        #self.domain.tight_slope_limiters = 1
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-
-        #Check contents
-        #Get NetCDF
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x']
-        y = fid.variables['y']
-        z = fid.variables['elevation']
-        time = fid.variables['time']
-        stage = fid.variables['stage']
-
-        #Check values
-        Q = self.domain.quantities['stage']
-        Q0 = Q.vertex_values[:,0]
-        Q1 = Q.vertex_values[:,1]
-        Q2 = Q.vertex_values[:,2]
-
-        A = stage[1,:]
-        assert num.allclose(A[0], min(Q2[0], Q1[1]))
-        assert num.allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
-        assert num.allclose(A[2], Q0[3])
-        assert num.allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
-
-        #Center point
-        assert num.allclose(A[4], min(Q1[0], Q2[1], Q0[2],
-                                      Q0[5], Q2[6], Q1[7]))
 
 
@@ -696 +639 @@
 
 
-    def Not_a_test_sww_DSG(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
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-
-        #Check contents
-        #Get NetCDF
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x']
-        y = fid.variables['y']
-        z = fid.variables['elevation']
-
-        volumes = fid.variables['volumes']
-        time = fid.variables['time']
-
-        # 2D
-        stage = fid.variables['stage']
-
-        X = x[:]
-        Y = y[:]
-        Z = z[:]
-        V = volumes[:]
-        T = time[:]
-        S = stage[:,:]
-
-#         print "****************************"
-#         print "X ",X
-#         print "****************************"
-#         print "Y ",Y
-#         print "****************************"
-#         print "Z ",Z
-#         print "****************************"
-#         print "V ",V
-#         print "****************************"
-#         print "Time ",T
-#         print "****************************"
-#         print "Stage ",S
-#         print "****************************"
-
-
-        fid.close()
-
-        #Cleanup
-        os.remove(sww.filename)
-
-
-
-    def test_export_grid(self):
-        """
-        test_export_grid(self):
-        Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        try:
-            os.remove('teg*.sww')
-        except:
-            pass
-
-        #Setup
-        self.domain.set_name('teg')
-
-        prjfile = self.domain.get_name() + '_elevation.prj'
-        ascfile = self.domain.get_name() + '_elevation.asc'
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 1.0)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-        fid.close()
-
-        #Export to ascii/prj files
-        sww2dem_batch(self.domain.get_name(),
-                quantities = 'elevation',
-                cellsize = cellsize,
-                verbose = self.verbose,
-                format = 'asc')
-
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), -i*cellsize - y)
-                
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-    def test_export_gridII(self):
-        """
-        test_export_gridII(self):
-        Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        try:
-            os.remove('teg*.sww')
-        except:
-            pass
-
-        #Setup
-        self.domain.set_name('tegII')
-
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 1.0)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep()
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep()
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]        
-
-        #print 'stage', stage
-        #print 'xmom', xmomentum
-        #print 'ymom', ymomentum        
-
-        fid.close()
-
-        #Export to ascii/prj files
-        if True:
-            sww2dem_batch(self.domain.get_name(),
-                        quantities = ['elevation', 'depth'],
-                        cellsize = cellsize,
-                        verbose = self.verbose,
-                        format = 'asc')
-
-        else:
-            sww2dem_batch(self.domain.get_name(),
-                quantities = ['depth'],
-                cellsize = cellsize,
-                verbose = self.verbose,
-                format = 'asc')
-
-
-            export_grid(self.domain.get_name(),
-                quantities = ['elevation'],
-                cellsize = cellsize,
-                verbose = self.verbose,
-                format = 'asc')
-
-        prjfile = self.domain.get_name() + '_elevation.prj'
-        ascfile = self.domain.get_name() + '_elevation.asc'
-        
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        #print "ascfile", ascfile
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                #print " -i*cellsize - y",  -i*cellsize - y
-                #print "float(L[i])", float(L[i])
-                assert num.allclose(float(L[i]), -i*cellsize - y)
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        
-        #Check asc file
-        ascfile = self.domain.get_name() + '_depth.asc'
-        prjfile = self.domain.get_name() + '_depth.prj'
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                #print " -i*cellsize - y",  -i*cellsize - y
-                #print "float(L[i])", float(L[i])                
-                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-
-    def test_export_gridIII(self):
-        """
-        test_export_gridIII
-        Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        try:
-            os.remove('teg*.sww')
-        except:
-            pass
-
-        #Setup
-        
-        self.domain.set_name('tegIII')
-
-        swwfile = self.domain.get_name() + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-        self.domain.set_quantity('stage', 1.0)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-        
-        sww = SWW_file(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep() #'stage')
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep() #'stage')
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, netcdf_mode_r)
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-        fid.close()
-
-        #Export to ascii/prj files
-        extra_name_out = 'yeah'
-        if True:
-            sww2dem_batch(self.domain.get_name(),
-                        quantities = ['elevation', 'depth'],
-                        extra_name_out = extra_name_out,
-                        cellsize = cellsize,
-                        verbose = self.verbose,
-                        format = 'asc')
-
-        else:
-            sww2dem_batch(self.domain.get_name(),
-                quantities = ['depth'],
-                cellsize = cellsize,
-                verbose = self.verbose,
-                format = 'asc')
-
-
-            sww2dem_batch(self.domain.get_name(),
-                quantities = ['elevation'],
-                cellsize = cellsize,
-                verbose = self.verbose,
-                format = 'asc')
-
-        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
-        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
-        
-        #Check asc file
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        #print "ascfile", ascfile
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                #print " -i*cellsize - y",  -i*cellsize - y
-                #print "float(L[i])", float(L[i])
-                assert num.allclose(float(L[i]), -i*cellsize - y)
-                
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        
-        #Check asc file
-        ascfile = self.domain.get_name() + '_depth_yeah.asc'
-        prjfile = self.domain.get_name() + '_depth_yeah.prj'
-        ascid = open(ascfile)
-        lines = ascid.readlines()
-        ascid.close()
-
-        L = lines[2].strip().split()
-        assert L[0].strip().lower() == 'xllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 308500)
-
-        L = lines[3].strip().split()
-        assert L[0].strip().lower() == 'yllcorner'
-        assert num.allclose(float(L[1].strip().lower()), 6189000)
-
-        #Check grid values
-        for j in range(5):
-            L = lines[6+j].strip().split()
-            y = (4-j) * cellsize
-            for i in range(5):
-                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
-
-        #Cleanup
-        os.remove(prjfile)
-        os.remove(ascfile)
-        os.remove(swwfile)
-
-    def test_export_grid_bad(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        try:
-            sww2dem_batch('a_small_round-egg',
-                        quantities = ['elevation', 'depth'],
-                        cellsize = 99,
-                        verbose = self.verbose,
-                        format = 'asc')
-        except IOError:
-            pass
-        else:
-            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
-
-
     def test_file_boundary_stsIV_sinewave_ordering(self):
         """test_file_boundary_stsIV_sinewave_ordering(self):
@@ -1195 +712 @@
         xmomentum = fid.variables['xmomentum'][:]
         ymomentum = fid.variables['ymomentum'][:]
-        elevation = fid.variables['elevation'][:]
-
-        #print stage
-        #print xmomentum
-        #print ymomentum
-        #print elevation
-        
-        
+        elevation = fid.variables['elevation'][:]       
 
         # Create beginnings of boundary polygon based on sts_boundary
@@ -1278 +788 @@
                                                    skip_initial_step=False)):
             temp_time[i]=domain_time.quantities['stage'].centroid_values[2]
-
-
-
-        #print temp_fbound
-        #print temp_time
-
-        #print domain_fbound.quantities['stage'].vertex_values
-        #print domain_time.quantities['stage'].vertex_values
         
         assert num.allclose(temp_fbound, temp_time)                
@@ -1306 +808 @@
         os.remove(meshname)
         
-        
-
-        
-        
-    def test_file_boundary_sts_time_limit(self):
-        """test_file_boundary_stsIV_sinewave_ordering(self):
-        Read correct points from ordering file and apply sts to boundary
-        This one uses a sine wave and compares to time boundary
-        
-        This one tests that times used can be limited by upper limit
-        """
-        
-        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
-        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
-        tide = 0.35
-        time_step_count = 50
-        time_step = 0.1
-        times_ref = num.arange(0, time_step_count*time_step, time_step)
-        
-        n=len(lat_long_points)
-        first_tstep=num.ones(n,num.int)
-        last_tstep=(time_step_count)*num.ones(n,num.int)
-        
-        gauge_depth=20*num.ones(n,num.float)
-        
-        ha1=num.ones((n,time_step_count),num.float)
-        ua1=3.*num.ones((n,time_step_count),num.float)
-        va1=2.*num.ones((n,time_step_count),num.float)
-        for i in range(n):
-            ha1[i]=num.sin(times_ref)
-        
-        
-        base_name, files = self.write_mux2(lat_long_points,
-                                           time_step_count, time_step,
-                                           first_tstep, last_tstep,
-                                           depth=gauge_depth,
-                                           ha=ha1,
-                                           ua=ua1,
-                                           va=va1)
-
-        # Write order file
-        file_handle, order_base_name = tempfile.mkstemp("")
-        os.close(file_handle)
-        os.remove(order_base_name)
-        d=","
-        order_file=order_base_name+'order.txt'
-        fid=open(order_file,'w')
-        
-        # Write Header
-        header='index, longitude, latitude\n'
-        fid.write(header)
-        indices=[3,0,1]
-        for i in indices:
-            line=str(i)+d+str(lat_long_points[i][1])+d+\
-                str(lat_long_points[i][0])+"\n"
-            fid.write(line)
-        fid.close()
-
-        sts_file=base_name
-        urs2sts(base_name, basename_out=sts_file,
-                ordering_filename=order_file,
-                mean_stage=tide,
-                verbose=False)
-        self.delete_mux(files)
-        
-        
-        
-        # Now read the sts file and check that values have been stored correctly.
-        fid = NetCDFFile(sts_file + '.sts')
-
-        # Check the time vector
-        times = fid.variables['time'][:]
-        starttime = fid.starttime[0]
-        #print times
-        #print starttime
-
-        # Check sts quantities
-        stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-        elevation = fid.variables['elevation'][:]
-
-        
-
-        # Create beginnings of boundary polygon based on sts_boundary
-        boundary_polygon = create_sts_boundary(base_name)
-        
-        os.remove(order_file)
-
-        # Append the remaining part of the boundary polygon to be defined by
-        # the user
-        bounding_polygon_utm=[]
-        for point in bounding_polygon:
-            zone,easting,northing=redfearn(point[0],point[1])
-            bounding_polygon_utm.append([easting,northing])
-
-        boundary_polygon.append(bounding_polygon_utm[3])
-        boundary_polygon.append(bounding_polygon_utm[4])
-
-        #print 'boundary_polygon', boundary_polygon
-        
-
-        assert num.allclose(bounding_polygon_utm,boundary_polygon)
-
-
-        extent_res=1000000
-        meshname = 'urs_test_mesh' + '.tsh'
-        interior_regions=None
-        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
-        
-        # have to change boundary tags from last example because now bounding
-        # polygon starts in different place.
-        create_mesh_from_regions(boundary_polygon,
-                                 boundary_tags=boundary_tags,
-                                 maximum_triangle_area=extent_res,
-                                 filename=meshname,
-                                 interior_regions=interior_regions,
-                                 verbose=False)
-        
-        domain_fbound = Domain(meshname)
-        domain_fbound.set_quantity('stage', tide)
-        
-        
-        Bf = File_boundary(sts_file+'.sts', 
-                           domain_fbound, 
-                           boundary_polygon=boundary_polygon)
-        time_vec = Bf.F.get_time()
-        assert num.allclose(Bf.F.starttime, starttime)
-        assert num.allclose(time_vec, times_ref)                                   
-        
-        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
-            Bf = File_boundary(sts_file+'.sts', 
-                               domain_fbound, 
-                               time_limit=time_limit+starttime,
-                               boundary_polygon=boundary_polygon)
-        
-            time_vec = Bf.F.get_time()
-            assert num.allclose(Bf.F.starttime, starttime)            
-            assert num.alltrue(time_vec < time_limit)
-            
-            
-        try:    
-            Bf = File_boundary(sts_file+'.sts', 
-                               domain_fbound, 
-                               time_limit=-1+starttime,
-                               boundary_polygon=boundary_polygon)            
-            time_vec = Bf.F.get_time()    
-            print time_vec    
-        except AssertionError:
-            pass
-        else:
-            raise Exception, 'Should have raised Exception here'
-
-    #### END TESTS FOR URS 2 SWW  ###
-
-
-    def test_triangulation(self):
-        # 
-        #  
-        
-        filename = tempfile.mktemp("_data_manager.sww")
-        outfile = NetCDFFile(filename, netcdf_mode_w)
-        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
-        volumes = (0,1,2)
-        elevation = [0,1,2]
-        new_origin = None
-        new_origin = Geo_reference(56, 0, 0)
-        times = [0, 10]
-        number_of_volumes = len(volumes)
-        number_of_points = len(points_utm)
-        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
-        sww.store_header(outfile, times, number_of_volumes,
-                         number_of_points, description='fully sick testing',
-                         verbose=self.verbose,sww_precision=netcdf_float)
-        sww.store_triangulation(outfile, points_utm, volumes,
-                                elevation,  new_origin=new_origin,
-                                verbose=self.verbose)       
-        outfile.close()
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        fid.close()
-
-        assert num.allclose(num.array(map(None, x,y)), points_utm)
-        os.remove(filename)
-
-        
-    def test_triangulationII(self):
-        # 
-        #  
-        
-        filename = tempfile.mktemp("_data_manager.sww")
-        outfile = NetCDFFile(filename, netcdf_mode_w)
-        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
-        volumes = (0,1,2)
-        elevation = [0,1,2]
-        new_origin = None
-        #new_origin = Geo_reference(56, 0, 0)
-        times = [0, 10]
-        number_of_volumes = len(volumes)
-        number_of_points = len(points_utm)
-        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
-        sww.store_header(outfile, times, number_of_volumes,
-                         number_of_points, description='fully sick testing',
-                         verbose=self.verbose,sww_precision=netcdf_float)
-        sww.store_triangulation(outfile, points_utm, volumes,
-                                new_origin=new_origin,
-                                verbose=self.verbose)
-        sww.store_static_quantities(outfile, elevation=elevation)                                
-                                
-        outfile.close()
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        results_georef = Geo_reference()
-        results_georef.read_NetCDF(fid)
-        assert results_georef == Geo_reference(DEFAULT_ZONE, 0, 0)
-        fid.close()
-
-        assert num.allclose(num.array(map(None, x,y)), points_utm)
-        os.remove(filename)
-
-        
-    def test_triangulation_new_origin(self):
-        # 
-        #  
-        
-        filename = tempfile.mktemp("_data_manager.sww")
-        outfile = NetCDFFile(filename, netcdf_mode_w)
-        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
-        volumes = (0,1,2)
-        elevation = [0,1,2]
-        new_origin = None
-        new_origin = Geo_reference(56, 1, 554354)
-        points_utm = new_origin.change_points_geo_ref(points_utm)
-        times = [0, 10]
-        number_of_volumes = len(volumes)
-        number_of_points = len(points_utm)
-        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
-        sww.store_header(outfile, times, number_of_volumes,
-                         number_of_points, description='fully sick testing',
-                         verbose=self.verbose,sww_precision=netcdf_float)
-        sww.store_triangulation(outfile, points_utm, volumes,
-                                elevation,  new_origin=new_origin,
-                                verbose=self.verbose)
-        outfile.close()
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        results_georef = Geo_reference()
-        results_georef.read_NetCDF(fid)
-        assert results_georef == new_origin
-        fid.close()
-
-        absolute = Geo_reference(56, 0,0)
-        assert num.allclose(num.array(
-            absolute.change_points_geo_ref(map(None, x,y),
-                                           new_origin)),points_utm)
-        
-        os.remove(filename)
-        
-    def test_triangulation_points_georeference(self):
-        # 
-        #  
-        
-        filename = tempfile.mktemp("_data_manager.sww")
-        outfile = NetCDFFile(filename, netcdf_mode_w)
-        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
-        volumes = (0,1,2)
-        elevation = [0,1,2]
-        new_origin = None
-        points_georeference = Geo_reference(56, 1, 554354)
-        points_utm = points_georeference.change_points_geo_ref(points_utm)
-        times = [0, 10]
-        number_of_volumes = len(volumes)
-        number_of_points = len(points_utm)
-        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
-        sww.store_header(outfile, times, number_of_volumes,
-                         number_of_points, description='fully sick testing',
-                         verbose=self.verbose,sww_precision=netcdf_float)
-        sww.store_triangulation(outfile, points_utm, volumes,
-                                elevation,  new_origin=new_origin,
-                                points_georeference=points_georeference,
-                                verbose=self.verbose)       
-        outfile.close()
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        results_georef = Geo_reference()
-        results_georef.read_NetCDF(fid)
-        assert results_georef == points_georeference
-        fid.close()
-
-        assert num.allclose(num.array(map(None, x,y)), points_utm)
-        os.remove(filename)
-        
-    def test_triangulation_2_geo_refs(self):
-        # 
-        #  
-        
-        filename = tempfile.mktemp("_data_manager.sww")
-        outfile = NetCDFFile(filename, netcdf_mode_w)
-        points_utm = num.array([[0.,0.],[1.,1.], [0.,1.]])
-        volumes = (0,1,2)
-        elevation = [0,1,2]
-        new_origin = Geo_reference(56, 1, 1)
-        points_georeference = Geo_reference(56, 0, 0)
-        points_utm = points_georeference.change_points_geo_ref(points_utm)
-        times = [0, 10]
-        number_of_volumes = len(volumes)
-        number_of_points = len(points_utm)
-        sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])        
-        sww.store_header(outfile, times, number_of_volumes,
-                         number_of_points, description='fully sick testing',
-                         verbose=self.verbose,sww_precision=netcdf_float)
-        sww.store_triangulation(outfile, points_utm, volumes,
-                                elevation,  new_origin=new_origin,
-                                points_georeference=points_georeference,
-                                verbose=self.verbose)       
-        outfile.close()
-        fid = NetCDFFile(filename)
-
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        results_georef = Geo_reference()
-        results_georef.read_NetCDF(fid)
-        assert results_georef == new_origin
-        fid.close()
-
-
-        absolute = Geo_reference(56, 0,0)
-        assert num.allclose(num.array(
-            absolute.change_points_geo_ref(map(None, x,y),
-                                           new_origin)),points_utm)
-        os.remove(filename)
-
  
     def test_points2polygon(self):  
@@ -1664 +826 @@
         assert (polygon == [[0.0, 0.0],[1.0, 0.0],[0.0, 1.0]])
         
-    
-    def test_csv2polygons(self):
-        """test_csv2polygons
-        """
-        
-        path = get_pathname_from_package('anuga.shallow_water')                
-        testfile = os.path.join(path, 'polygon_values_example.csv')                
-
-        polygons, values = load_csv_as_polygons(testfile, 
-                                        value_name='floors')
-
-        assert len(polygons) == 7, 'Must have seven polygons'
-        assert len(values) == 7, 'Must have seven values'
-            
-        # Known floor values
-        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
-        
-        # Known polygon values
-        known_polys = {'1': [[422681.61,871117.55],
-                             [422691.02,871117.60],
-                             [422690.87,871084.23],
-                             [422649.36,871081.85],
-                             [422649.36,871080.39],
-                             [422631.86,871079.50],
-                             [422631.72,871086.75],
-                             [422636.75,871087.20],
-                             [422636.75,871091.50],
-                             [422649.66,871092.09],
-                             [422649.83,871084.91],
-                             [422652.94,871084.90],
-                             [422652.84,871092.39],
-                             [422681.83,871093.73],
-                             [422681.61,871117.55]],
-                       '2': [[422664.22,870785.46],
-                             [422672.48,870780.14],
-                             [422668.17,870772.62],
-                             [422660.35,870777.17],
-                             [422664.22,870785.46]],
-                       '3': [[422661.30,871215.06],
-                             [422667.50,871215.70],
-                             [422668.30,871204.86],
-                             [422662.21,871204.33],
-                             [422661.30,871215.06]],
-                       '4': [[422473.44,871191.22],
-                             [422478.33,871192.26],
-                             [422479.52,871186.03],
-                             [422474.78,871185.14],
-                             [422473.44,871191.22]],
-                       '5': [[422369.69,871049.29],
-                             [422378.63,871053.58],
-                             [422383.91,871044.51],
-                             [422374.97,871040.32],
-                             [422369.69,871049.29]],
-                       '8': [[422730.56,871203.13],
-                             [422734.10,871204.90],
-                             [422735.26,871202.18],
-                             [422731.87,871200.58],
-                             [422730.56,871203.13]],
-                       '9': [[422659.85,871213.80],
-                             [422660.91,871210.97],
-                             [422655.42,871208.85],
-                             [422654.36,871211.68],
-                             [422659.85,871213.80]]
-                       }        
-        
-
-        
-                
-        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
-            assert id in polygons.keys()
-            assert id in values.keys()            
-
-            assert int(values[id]) == int(floors[id])
-            assert len(polygons[id]) == len(known_polys[id])
-            assert num.allclose(polygons[id], known_polys[id])
-
-
-    def test_csv2polygons_with_clipping(self):
-        """test_csv2polygons with optional clipping
-        """
-        #FIXME(Ole): Not Done!!
-        
-        path = get_pathname_from_package('anuga.shallow_water')                
-        testfile = os.path.join(path, 'polygon_values_example.csv')                
-
-        polygons, values = load_csv_as_polygons(testfile, 
-                                        value_name='floors',
-                                        clipping_polygons=None)
-
-        assert len(polygons) == 7, 'Must have seven polygons'
-        assert len(values) == 7, 'Must have seven values'
-            
-        # Known floor values
-        floors = {'1': 2, '2': 0, '3': 1, '4': 2, '5': 0, '8': 1, '9': 1}
-        
-        # Known polygon values
-        known_polys = {'1': [[422681.61,871117.55],
-                             [422691.02,871117.60],
-                             [422690.87,871084.23],
-                             [422649.36,871081.85],
-                             [422649.36,871080.39],
-                             [422631.86,871079.50],
-                             [422631.72,871086.75],
-                             [422636.75,871087.20],
-                             [422636.75,871091.50],
-                             [422649.66,871092.09],
-                             [422649.83,871084.91],
-                             [422652.94,871084.90],
-                             [422652.84,871092.39],
-                             [422681.83,871093.73],
-                             [422681.61,871117.55]],
-                       '2': [[422664.22,870785.46],
-                             [422672.48,870780.14],
-                             [422668.17,870772.62],
-                             [422660.35,870777.17],
-                             [422664.22,870785.46]],
-                       '3': [[422661.30,871215.06],
-                             [422667.50,871215.70],
-                             [422668.30,871204.86],
-                             [422662.21,871204.33],
-                             [422661.30,871215.06]],
-                       '4': [[422473.44,871191.22],
-                             [422478.33,871192.26],
-                             [422479.52,871186.03],
-                             [422474.78,871185.14],
-                             [422473.44,871191.22]],
-                       '5': [[422369.69,871049.29],
-                             [422378.63,871053.58],
-                             [422383.91,871044.51],
-                             [422374.97,871040.32],
-                             [422369.69,871049.29]],
-                       '8': [[422730.56,871203.13],
-                             [422734.10,871204.90],
-                             [422735.26,871202.18],
-                             [422731.87,871200.58],
-                             [422730.56,871203.13]],
-                       '9': [[422659.85,871213.80],
-                             [422660.91,871210.97],
-                             [422655.42,871208.85],
-                             [422654.36,871211.68],
-                             [422659.85,871213.80]]
-                       }        
-        
-
-        
-                
-        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
-            assert id in polygons.keys()
-            assert id in values.keys()            
-
-            assert int(values[id]) == int(floors[id])
-            assert len(polygons[id]) == len(known_polys[id])
-            assert num.allclose(polygons[id], known_polys[id])
-
-
-
-
-    
-    def test_csv2building_polygons(self):
-        """test_csv2building_polygons
-        """
-        
-        path = get_pathname_from_package('anuga.shallow_water')                
-        testfile = os.path.join(path, 'polygon_values_example.csv')                
-
-        polygons, values = load_csv_as_building_polygons(testfile, 
-                                                 floor_height=3)
-
-        assert len(polygons) == 7, 'Must have seven polygons'
-        assert len(values) == 7, 'Must have seven values'
-            
-        # Known floor values
-        floors = {'1': 6, '2': 0, '3': 3, '4': 6, '5': 0, '8': 3, '9': 3}
-        
-                
-        for id in ['1', '2', '3', '4', '5' ,'8' ,'9']:
-            assert id in polygons.keys()
-            assert id in values.keys()            
-
-            assert float(values[id]) == float(floors[id])
 
 

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py

@@ -6 +6 @@
 import tempfile
 
+import anuga
+
 from anuga.config import g, epsilon
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
@@ -2085 +2087 @@
         can be controlled this way.
         """
-
-        #---------------------------------------------------------------------
-        # Import necessary modules
-        #---------------------------------------------------------------------
-        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water import Reflective_boundary
-        from anuga.shallow_water import Dirichlet_boundary
-        from anuga.shallow_water import Time_boundary
-
+        
         #---------------------------------------------------------------------
         # Setup computational domain
@@ -2108 +2101 @@
                                                        len1=length, 
                                                        len2=width)
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.set_name('channel_variable_test')  # Output name
         domain.set_quantities_to_be_stored({'elevation': 2,

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py

@@ -9 +9 @@
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_conservation.py

@@ -6 +6 @@
 import tempfile
 
+import anuga
+
 from anuga.config import g, epsilon
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
@@ -60 +62 @@
 
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -88 +90 @@
         """
 
-        from mesh_factory import rectangular
-
         # Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = anuga.rectangular(6, 6)
 
         # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.smooth = False
         domain.default_order = 2
@@ -107 +107 @@
 
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -140 +140 @@
 
         # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.smooth = False
 
@@ -163 +163 @@
 
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -200 +200 @@
         """
 
-        from mesh_factory import rectangular
-
         # Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = anuga.rectangular(6, 6)
 
         # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.smooth = False
         domain.default_order = 2
@@ -232 +230 @@
 
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -271 +269 @@
         """
 
-        from mesh_factory import rectangular
-
         # Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = anuga.rectangular(6, 6)
 
         # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.smooth = False
         domain.default_order = 2
@@ -290 +286 @@
 
         # Boundary conditions (reflective everywhere)
-        Br = Reflective_boundary(domain)
-        Bleft = Dirichlet_boundary([0.5, 0, 0])
-        Bright = Dirichlet_boundary([0.1, 0, 0])
+        Br = anuga.Reflective_boundary(domain)
+        Bleft = anuga.Dirichlet_boundary([0.5, 0, 0])
+        Bright = anuga.Dirichlet_boundary([0.1, 0, 0])
         domain.set_boundary({'left': Bleft, 'right': Bright,
                              'top': Br, 'bottom': Br})
@@ -362 +358 @@
 
         # Boundaries
-        R = Reflective_boundary(domain)
+        R = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': R, 'right': R, 'top':R, 'bottom': R})
 
@@ -384 +380 @@
         Test that total volume can be computed correctly
         """            
-
-        #----------------------------------------------------------------------
-        # Import necessary modules
-        #----------------------------------------------------------------------
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from anuga.shallow_water import Domain
 
         #----------------------------------------------------------------------
@@ -449 +438 @@
         verbose = False
 
-        #----------------------------------------------------------------------
-        # Import necessary modules
-        #----------------------------------------------------------------------
-
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.forcing import Inflow
-        from anuga.shallow_water.data_manager \
-                import get_flow_through_cross_section
+
 
         #----------------------------------------------------------------------
@@ -500 +478 @@
         #----------------------------------------------------------------------
 
-        Br = Reflective_boundary(domain)      # Solid reflective wall
+        Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
                 
         # Constant flow in and out of domain
         # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 
-        Bi = Dirichlet_boundary([d, uh, vh]) 
-        Bo = Dirichlet_boundary([d, uh, vh])
+        Bi = anuga.Dirichlet_boundary([d, uh, vh]) 
+        Bo = anuga.Dirichlet_boundary([d, uh, vh])
 
         domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
@@ -544 +522 @@
         
 
-        #---------------------------------------------------------------------
-        # Import necessary modules
-        #---------------------------------------------------------------------
-        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.forcing import Inflow
-        from anuga.shallow_water.data_manager import get_flow_through_cross_section
-
         #----------------------------------------------------------------------
         # Setup computational domain
@@ -569 +537 @@
 
 
-        domain = Domain(points, vertices, boundary)   
+        domain = anuga.Domain(points, vertices, boundary)   
         domain.set_name('Inflow_volume_test')              # Output name
                 
@@ -593 +561 @@
 
         # Fixed Flowrate onto Area 
-        fixed_inflow = Inflow(domain,
+        fixed_inflow = anuga.Inflow(domain,
                               center=(10.0, 10.0),
                               radius=5.00,
@@ -604 +572 @@
         #----------------------------------------------------------------------
 
-        Br = Reflective_boundary(domain) # Solid reflective wall
+        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -646 +614 @@
         
 
-        #---------------------------------------------------------------------
-        # Import necessary modules
-        #---------------------------------------------------------------------
-        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.shallow_water_domain import Rainfall
-        from anuga.shallow_water.data_manager import get_flow_through_cross_section
-
         #----------------------------------------------------------------------
         # Setup computational domain
@@ -666 +624 @@
         
 
-        points, vertices, boundary = rectangular_cross(int(length/dx), 
+        points, vertices, boundary = anuga.rectangular_cross(int(length/dx), 
                                                        int(width/dy),
                                                        len1=length, len2=width)
 
 
-        domain = Domain(points, vertices, boundary)   
+        domain = anuga.Domain(points, vertices, boundary)   
         domain.set_name('Rain_volume_test')              # Output name
                 
@@ -695 +653 @@
 
         # Fixed rain onto small circular area 
-        fixed_rain = Rainfall(domain,
+        fixed_rain = anuga.Rainfall(domain,
                               center=(10.0, 10.0),
                               radius=5.00,
@@ -706 +664 @@
         #----------------------------------------------------------------------
 
-        Br = Reflective_boundary(domain) # Solid reflective wall
+        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -755 +713 @@
         
 
-        #---------------------------------------------------------------------
-        # Import necessary modules
-        #---------------------------------------------------------------------
-        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.shallow_water_domain import Rainfall
-        from anuga.shallow_water.data_manager import get_flow_through_cross_section
-
         #----------------------------------------------------------------------
         # Setup computational domain
@@ -815 +763 @@
         #----------------------------------------------------------------------
 
-        Br = Reflective_boundary(domain) # Solid reflective wall
-        Bt = Transmissive_stage_zero_momentum_boundary(domain)
-        Bd = Dirichlet_boundary([-10, 0, 0])
+        Br = anuga.Reflective_boundary(domain) # Solid reflective wall
+        Bt = anuga.Transmissive_stage_zero_momentum_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([-10, 0, 0])
         domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
 

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_cross_sections.py

@@ -6 +6 @@
 import tempfile
 
+import anuga
+
 from anuga.config import g, epsilon
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
-from anuga.abstract_2d_finite_volumes.quantity import Quantity
-from anuga.geospatial_data.geospatial_data import Geospatial_data
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 
 from anuga.utilities.system_tools import get_pathname_from_package
@@ -76 +75 @@
         uh = u*h
 
-        Br = Reflective_boundary(domain)     # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
+        Br = anuga.Reflective_boundary(domain)     # Side walls
+        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet: 
 
 
@@ -161 +160 @@
         uh = u*h
 
-        Br = Reflective_boundary(domain)       # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
+        Br = anuga.Reflective_boundary(domain)     # Side walls
+        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
 
         # Initial conditions
@@ -252 +251 @@
         uh = u*h
 
-        Br = Reflective_boundary(domain)       # Side walls
-        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
+        Br = anuga.Reflective_boundary(domain)     # Side walls
+        Bd = anuga.Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
 
         # Initial conditions

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py

@@ -9 +9 @@
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

@@ -5 +5 @@
 from math import pi, sqrt
 import tempfile
+import anuga
 
 from anuga.config import g, epsilon
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
@@ -394 +395 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        domain = Domain(points, vertices)
+        domain = anuga.Domain(points, vertices)
 
         #Flat surface with 1m of water
@@ -401 +402 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -408 +409 @@
         phi = 135
         domain.forcing_terms = []
-        domain.forcing_terms.append(Wind_stress(s, phi))
+        domain.forcing_terms.append(anuga.Wind_stress(s, phi))
 
         domain.compute_forcing_terms()
@@ -450 +451 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -459 +460 @@
         phi = 135
         domain.forcing_terms = []
-        domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
+        domain.forcing_terms.append(anuga.Wind_stress(s=speed, phi=angle))
 
         domain.compute_forcing_terms()
@@ -522 +523 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -543 +544 @@
         fid.close()
 
-        # Convert ASCII file to NetCDF (Which is what we really like!)
-        from data_manager import timefile2netcdf
-
-        timefile2netcdf(filename)
+        anuga.timefile2netcdf(filename+'.txt', filename+'.tms')
         os.remove(filename + '.txt')
 
@@ -616 +614 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -623 +621 @@
         # Write wind stress file (ensure that domain.time is covered)
         # Take x=1 and y=0
-        filename = 'test_windstress_from_file'
+        filename = 'test_windstress_from_file.txt'
+        file_out = 'test_windstress_from_file.tms'
         start = time.mktime(time.strptime('2000', '%Y'))
-        fid = open(filename + '.txt', 'w')
+        fid = open(filename, 'w')
         dt = 0.5    # Half second interval
         t = 0.0
@@ -635 +634 @@
         fid.close()
 
-        # Convert ASCII file to NetCDF (Which is what we really like!)
-        from data_manager import timefile2netcdf
-
-        timefile2netcdf(filename, time_as_seconds=True)
-        os.remove(filename + '.txt')
+        anuga.timefile2netcdf(filename, file_out, time_as_seconds=True)
+        os.remove(filename)
 
         # Setup wind stress
-        F = file_function(filename + '.tms',
+        F = file_function(file_out,
                           quantities=['Attribute0', 'Attribute1'])
-        os.remove(filename + '.tms')
-
-        W = Wind_stress(F)
+        os.remove(file_out)
+
+        W = anuga.Wind_stress(F)
 
         domain.forcing_terms = []
@@ -709 +705 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -718 +714 @@
 
         try:
-            domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
+            domain.forcing_terms.append(anuga.Wind_stress(s=scalar_func_list,
                                                     phi=angle))
         except AssertionError:
@@ -763 +759 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
         # Setup only one forcing term, constant rainfall
         domain.forcing_terms = []
-        domain.forcing_terms.append(Rainfall(domain, rate=2.0))
+        domain.forcing_terms.append(anuga.Rainfall(domain, rate=2.0))
 
         domain.compute_forcing_terms()
@@ -795 +791 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -801 +797 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
+        R = anuga.Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
 
         assert num.allclose(R.exchange_area, 2)
@@ -826 +822 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        domain = Domain(points, vertices)
+        domain = anuga.Domain(points, vertices)
 
         # Flat surface with 1m of water
@@ -833 +829 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -839 +835 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain,
+        R = anuga.Rainfall(domain,
                      rate=lambda t: 3*t + 7,
                      polygon = [[1,1], [2,1], [2,2], [1,2]])
@@ -870 +866 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        domain = Domain(points, vertices)
+        domain = anuga.Domain(points, vertices)
 
         # Flat surface with 1m of water
@@ -877 +873 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -883 +879 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain,
+        R = anuga.Rainfall(domain,
                      rate=lambda t: 3*t + 7,
                      polygon=rainfall_poly)                     
@@ -943 +939 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -949 +945 @@
         # restricted to a polygon enclosing triangle #1 (bce)
         domain.forcing_terms = []
-        R = Rainfall(domain,
+        R = anuga.Rainfall(domain,
                      rate=lambda t: 3*t + 7,
                      polygon=rainfall_poly)
@@ -1000 +996 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -1015 +1011 @@
 
         domain.forcing_terms = []
-        R = Rainfall(domain,
+        R = anuga.Rainfall(domain,
                      rate=main_rate,
                      polygon = [[1,1], [2,1], [2,2], [1,2]],
@@ -1068 +1064 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -1083 +1079 @@
 
         domain.forcing_terms = []
-        R = Rainfall(domain,
+        R = anuga.Rainfall(domain,
                      rate=main_rate,
                      polygon=[[1,1], [2,1], [2,2], [1,2]],
@@ -1120 +1116 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -1127 +1123 @@
         domain.forcing_terms = []
         
-        I = Inflow(domain, rate=2.0, center=(1,1), radius=1)
+        I = anuga.Inflow(domain, rate=2.0, center=(1,1), radius=1)
         domain.forcing_terms.append(I)        
         domain.compute_forcing_terms()
@@ -1154 +1150 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        domain = Domain(points, vertices)
+        domain = anuga.Domain(points, vertices)
 
         # Flat surface with 1m of water
@@ -1161 +1157 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -1208 +1204 @@
         domain.set_quantity('friction', 0)
 
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'exterior': Br})
 
@@ -1352 +1348 @@
         verbose = False
         
-
-        #----------------------------------------------------------------------
-        # Import necessary modules
-        #----------------------------------------------------------------------
-
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.forcing import Inflow
-        from anuga.shallow_water.data_manager \
-                import get_flow_through_cross_section
-        from anuga.abstract_2d_finite_volumes.util \
-                import sww2csv_gauges, csv2timeseries_graphs
-
         #----------------------------------------------------------------------
         # Setup computational domain
@@ -1414 +1394 @@
 
                 # Fixed Flowrate onto Area 
-                fixed_inflow = Inflow(domain,
+                fixed_inflow = anuga.Inflow(domain,
                                       center=(10.0, 10.0),
                                       radius=5.00,
@@ -1423 +1403 @@
                     domain.forcing_terms.append(fixed_inflow)
                 
-                ref_flow = fixed_inflow.rate*number_of_inflows
+                ref_flow = fixed_inflow.rate*number_of_inflowsg
 
                 # Compute normal depth on plane using Mannings equation

trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_reporting.py

@@ -6 +6 @@
 import tempfile
 
+import anuga
+
 from anuga.config import g, epsilon
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
 from anuga.utilities.numerical_tools import mean
-from anuga.utilities.polygon import is_inside_polygon
+from anuga.geometry.polygon import is_inside_polygon
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
@@ -77 +79 @@
         # Setup boundary conditions
         #--------------------------------------------------------------
-        Br = Reflective_boundary(domain)              # Reflective wall
+        Br = anuga.Reflective_boundary(domain)       # Reflective wall
         # Constant inflow
-        Bd = Dirichlet_boundary([final_runup_height, 0, 0])
+        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])
 
         # All reflective to begin with (still water)