Changeset 1669

Timestamp:

Aug 2, 2005, 12:36:06 PM (20 years ago)

Author:

hitchman

Message:

added test in decimate_dem to deal sensibly, if roughly, with NODATA values

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (2 diffs)
• test_data_manager.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -2402 +2402 @@
     """
 
-    #FIXME: does not yet deal sensibly with NODATA values, they are
-    #       included in calculations
-
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, zeros, sum, reshape
+    from Numeric import Float, zeros, sum, reshape, equal
 
     root = basename_in
@@ -2493 +2490 @@
         for j in range(ncols_new):
             tcol = j * cellsize_ratio
-            telev[local_index] = \
-             sum(sum(dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil] * stencil))
+            tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil]
+
+            #if dem contains 1 or more NODATA_values set value in
+            #decimated dem to NODATA_value, else compute decimated
+            #value using stencil
+            if sum(sum(equal(tmp, NODATA_value))) > 0:
+                telev[local_index] = NODATA_value
+            else:
+                telev[local_index] = sum(sum(tmp * stencil))
+
             global_index += 1
             local_index += 1

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -1902 +1902 @@
         os.remove(root + '_100.dem')
 
+    def test_decimate_dem_NODATA(self):
+        """Test decimation of dem file that includes NODATA values
+        """
+
+        import os
+        from Numeric import ones, allclose, Float, arange, reshape
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test dem file
+        root = 'decdemtest'
+
+        filename = root + '.dem'
+        fid = NetCDFFile(filename, 'w')
+
+        fid.institution = 'Geoscience Australia'
+        fid.description = 'NetCDF DEM format for compact and portable ' +\
+                          'storage of spatial point data'
+
+        nrows = 15
+        ncols = 18
+        NODATA_value = -9999
+
+        fid.ncols = ncols
+        fid.nrows = nrows
+        fid.xllcorner = 2000.5
+        fid.yllcorner = 3000.5
+        fid.cellsize = 25
+        fid.NODATA_value = NODATA_value
+
+        fid.zone = 56
+        fid.false_easting = 0.0
+        fid.false_northing = 0.0
+        fid.projection = 'UTM'
+        fid.datum = 'WGS84'
+        fid.units = 'METERS'
+
+        fid.createDimension('number_of_points', nrows*ncols)
+
+        fid.createVariable('elevation', Float, ('number_of_points',))
+
+        elevation = fid.variables['elevation']
+
+        #generate initial elevation values
+        elevation_tmp = (arange(nrows*ncols))
+        #add some NODATA values
+        elevation_tmp[0]   = NODATA_value
+        elevation_tmp[95]  = NODATA_value
+        elevation_tmp[188] = NODATA_value
+        elevation_tmp[189] = NODATA_value
+        elevation_tmp[190] = NODATA_value
+        elevation_tmp[209] = NODATA_value
+        elevation_tmp[252] = NODATA_value
+
+        elevation[:] = elevation_tmp
+
+        fid.close()
+
+        #generate the elevation values expected in the decimated file
+        ref_elevation = [NODATA_value,
+                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
+                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
+                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
+                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
+                         NODATA_value,
+                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
+                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
+                         (144+145+146+162+163+164+180+181+182) / 9.0,
+                         (148+149+150+166+167+168+184+185+186) / 9.0,
+                         NODATA_value,
+                         (156+157+158+174+175+176+192+193+194) / 9.0,
+                         NODATA_value,
+                         (220+221+222+238+239+240+256+257+258) / 9.0,
+                         (224+225+226+242+243+244+260+261+262) / 9.0,
+                         (228+229+230+246+247+248+264+265+266) / 9.0]
+
+        #generate a stencil for computing the decimated values
+        stencil = ones((3,3), Float) / 9.0
+
+        decimate_dem(root, stencil=stencil, cellsize_new=100)
+
+        #Open decimated NetCDF file
+        fid = NetCDFFile(root + '_100.dem', 'r')
+
+        # Get decimated elevation
+        elevation = fid.variables['elevation']
+
+        #Check values
+        assert allclose(elevation, ref_elevation)
+
+        #Cleanup
+        fid.close()
+
+        os.remove(root + '.dem')
+        os.remove(root + '_100.dem')
+
 
 #-------------------------------------------------------------
@@ -1907 +2002 @@
     #suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
-    suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
+    suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
     runner = unittest.TextTestRunner()
     runner.run(suite)