Changeset 4462

Timestamp:

May 21, 2007, 10:15:58 AM (17 years ago)

Author:

duncan

Message:

adding export_grid

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (7 diffs)
• test_data_manager.py (modified) (6 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -86 +86 @@
      pmesh_to_domain_instance
 from anuga.abstract_2d_finite_volumes.util import get_revision_number
+
+# formula mappings
+
+quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
+                    'depth':'stage-elevation',
+                    'speed': \
+ '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
+
+
     
 def make_filename(s):
@@ -1597 +1606 @@
     geo.export_points_file(ptsname)
 
-
+def export_grid(basename_in, extra_name_out = None,
+                quantities = None, # defaults to elevation
+                timestep = None,
+                reduction = None,
+                cellsize = 10,
+                NODATA_value = -9999,
+                easting_min = None,
+                easting_max = None,
+                northing_min = None,
+                northing_max = None,
+                verbose = False,
+                origin = None,
+                datum = 'WGS84',
+                format = 'ers'):
+    """
+    
+    Wrapper for sww2dem. - see sww2dem to find out what most of the
+    parameters do.
+
+    Quantities is a list of quantities.  Each quantity will be
+    calculated for each sww file.
+
+    This returns the basenames of the files returned, which is made up
+    of the dir and all of the file name, except the extension.
+
+    This function returns the names of the files produced.
+    """
+    
+    if quantities is None:
+        quantities = ['elevation']
+        
+    if type(quantities) is str:
+            quantities = [quantities]
+
+    # How many sww files are there?
+    dir, base = os.path.split(basename_in)
+    dir_ls = os.listdir(dir)
+    interate_over = [x[:-4] for x in dir_ls if base in x and x[-4:] == '.sww']
+    #print "interate_over", interate_over
+    
+    files_out = []
+    for sww_file in interate_over:
+        for quantity in quantities:
+            if extra_name_out is None:
+                basename_out = sww_file + '_' + quantity
+            else:
+                basename_out = sww_file + '_' + quantity + '_' \
+                               + extra_name_out
+            #print "basename_out", basename_out
+        
+            file_out = sww2dem(sww_file, basename_out,
+                               quantity, 
+                               timestep,
+                               reduction,
+                               cellsize,
+                               NODATA_value,
+                               easting_min,
+                               easting_max,
+                               northing_min,
+                               northing_max,
+                               verbose,
+                               origin,
+                               datum,
+                               format)
+            files_out.append(file_out)
+    #print "basenames_out after",basenames_out 
+    return files_out
+    
 def sww2dem(basename_in, basename_out = None,
-            quantity = None,
+            quantity = None, # defaults to elevation
             timestep = None,
             reduction = None,
@@ -1611 +1687 @@
             origin = None,
             datum = 'WGS84',
-        format = 'ers'):
+            format = 'ers'):
 
     """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
@@ -1643 +1719 @@
     The parameter quantity must be the name of an existing quantity or
     an expression involving existing quantities. The default is
-    'elevation'.
+    'elevation'. Quantity is not a list of quantities.
 
     if timestep (an index) is given, output quantity at that timestep
@@ -1670 +1746 @@
     if quantity is None:
         quantity = 'elevation'
-
+        
     if reduction is None:
         reduction = max
@@ -1677 +1753 @@
         basename_out = basename_in + '_%s' %quantity
 
+    if quantity_formula.has_key(quantity):
+        quantity = quantity_formula[quantity]
+        
     swwfile = basename_in + '.sww'
     demfile = basename_out + '.' + format
@@ -1956 +2035 @@
         ascid.close()
         fid.close()
+
+        return basename_out
 
 #Backwards compatibility

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -1105 +1105 @@
 
 
-
-    def test_sww2dem_asc_elevation(self):
+    def test_sww2dem_asc_elevation_depth(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
@@ -1126 +1125 @@
         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)
@@ -1151 +1151 @@
         stage = fid.variables['stage'][:]
 
+        fid.close()
 
         #Export to ascii/prj files
@@ -1235 +1236 @@
             for i in range(5):
                 assert allclose(float(L[i]), -i*cellsize - y)
-
-
-        fid.close()
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.get_name(),
+                quantity = 'depth',
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+        
+        #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[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
+
 
         #Cleanup
@@ -1245 +1293 @@
 
 
+    def 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 Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        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.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 = get_dataobject(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, '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
+        export_grid(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 allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert 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 allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def 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 Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        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 = get_dataobject(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, '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
+        if True:
+            export_grid(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            export_grid(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 allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert 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 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 allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert 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 allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+    def test_export_gridIII(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.set_name('datatest')
+
+        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 = get_dataobject(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, '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:
+            export_grid(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        extra_name_out = extra_name_out,
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            export_grid(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_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 allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert 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 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 allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert 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 allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_export_grid_parallel(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        base_name = 'tegp'
+        #Setup
+        self.domain.set_name(base_name+'_P0_8')
+        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 = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+        self.domain.evolve_to_end(finaltime = 0.0001)
+        #Setup
+        self.domain.set_name(base_name+'_P1_8')
+        swwfile2 = self.domain.get_name() + '.sww'
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+        self.domain.evolve_to_end(finaltime = 0.0002)
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, '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'
+        export_grid(base_name,
+                    quantities = ['elevation', 'depth'],
+                    extra_name_out = extra_name_out,
+                    cellsize = cellsize,
+                    verbose = self.verbose,
+                    format = 'asc')
+
+        prjfile = base_name + '_P0_8_elevation_yeah.prj'
+        ascfile = base_name + '_P0_8_elevation_yeah.asc'       
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #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 allclose(float(L[i]), -i*cellsize - y)               
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        prjfile = base_name + '_P1_8_elevation_yeah.prj'
+        ascfile = base_name + '_P1_8_elevation_yeah.asc'       
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #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 allclose(float(L[i]), -i*cellsize - y)               
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+        #Check asc file
+        ascfile = base_name + '_P0_8_depth_yeah.asc'
+        prjfile = base_name + '_P0_8_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+
+        #Check asc file
+        ascfile = base_name + '_P1_8_depth_yeah.asc'
+        prjfile = base_name + '_P1_8_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert allclose(float(L[i]), 1 - (-i*cellsize - y))
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile2)
 
     def test_sww2dem_larger(self):
@@ -6187 +6687 @@
     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sww_origin')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_header')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww_range')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid_parallel')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_grid')
     if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
         Test_Data_Manager.verbose=True