Changeset 1867

Timestamp:

Oct 5, 2005, 5:37:19 PM (19 years ago)

Author:

ole

Message:

Added new unit test for bounding box in sww2dem - this functionality has not been implemented yet and the test is still disabled.

File:

• inundation/pyvolution/test_data_manager.py (modified) (2 diffs)

inundation/pyvolution/test_data_manager.py

@@ -17 +17 @@
         import time
         from mesh_factory import rectangular
-
 
         #Create basic mesh
@@ -875 +874 @@
         os.remove(ascfile)
         os.remove(swwfile)
+
+
+
+    def test_sww2dem_larger(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView. Here:
+
+        ncols         11
+        nrows         11
+        xllcorner     308500
+        yllcorner     6189000
+        cellsize      10.000000
+        NODATA_value  -9999
+        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200 
+         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 
+         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 
+         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 
+         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 
+         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 
+         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
+         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
+         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
+         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
+           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
+        
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+
+        from mesh_factory import rectangular
+
+        #Create basic mesh (100m x 100m)
+        points, vertices, boundary = rectangular(2, 2, 100, 100)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        domain.filename = 'datatest'
+
+        prjfile = domain.filename + '_elevation.prj'
+        ascfile = domain.filename + '_elevation.asc'
+        swwfile = domain.filename + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        #
+        domain.set_quantity('elevation', lambda x,y: -x-y)
+        domain.set_quantity('stage', 0)        
+
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        #
+        sww = get_dataobject(domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 10  #10m grid
+
+        
+        #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'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(domain.filename,
+                quantity = 'elevation',
+                cellsize = cellsize,
+                verbose = False,
+                format = 'asc')
+        #dem2pts(root, easting_min=2010.0, easting_max=2110.0,
+        #        northing_min=3035.0, northing_max=3125.5)
+        
+        
+                
+        #Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '11'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '11'
+
+        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 (FIXME: Use same strategy for other sww2dem tests)
+        for i, line in enumerate(lines[6:]):
+            for j, value in enumerate( line.split() ):
+                #assert allclose(float(value), -(10-i+j)*cellsize)
+                assert float(value) == -(10-i+j)*cellsize
+        
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+
+    #FIXME: TODO
+    def xtest_sww2dem_boundingbox(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView.
+        This will test that mesh can be restricted by bounding box
+
+        Original extent is 100m x 100m:
+
+        Eastings:   308500 -  308600 
+        Northings: 6189000 - 6189100
+
+        Bounding box changes this to the 50m x 50m square defined by
+
+        Eastings:   308530 -  308570 
+        Northings: 6189050 - 6189100
+
+        The cropped values should be
+
+         -130 -140 -150 -160 -170  
+         -120 -130 -140 -150 -160  
+         -110 -120 -130 -140 -150  
+         -100 -110 -120 -130 -140  
+          -90 -100 -110 -120 -130  
+          -80  -90 -100 -110 -120
+
+        and the new lower reference point should be   
+        Eastings:   308530 
+        Northings: 6189050
+        
+        Original dataset is the same as in test_sww2dem_larger()
+        
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+
+        from mesh_factory import rectangular
+
+        #Create basic mesh (100m x 100m)
+        points, vertices, boundary = rectangular(2, 2, 100, 100)
+
+        #Create shallow water domain
+        domain = Domain(points, vertices, boundary)
+        domain.default_order = 2
+
+        domain.filename = 'datatest'
+
+        prjfile = domain.filename + '_elevation.prj'
+        ascfile = domain.filename + '_elevation.asc'
+        swwfile = domain.filename + '.sww'
+
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.geo_reference = Geo_reference(56, 308500, 6189000)
+
+        #
+        domain.set_quantity('elevation', lambda x,y: -x-y)
+        domain.set_quantity('stage', 0)        
+
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+
+        #
+        sww = get_dataobject(domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 10  #10m grid
+
+        
+        #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'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(domain.filename,
+                quantity = 'elevation',
+                cellsize = cellsize,
+                easting_min = 308530,
+                easting_max = 308570,
+                northing_min = 6189050, 
+                northing_max = 6189100,                                
+                verbose = False,
+                format = 'asc')
+        
+                
+        #Check prj (meta data)
+        prjid = open(prjfile)
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+
+
+        #Check asc file
+        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() == '6'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308530)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189050)
+
+        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 i, line in enumerate(lines[6:]):
+            for j, value in enumerate( line.split() ):
+                #assert float(value) == -(10-i+j)*cellsize                
+                assert float(value) == -(10-i+j+3)*cellsize
+        
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+