Changeset 1666

Timestamp:

Aug 1, 2005, 12:57:53 PM (19 years ago)

Author:

hitchman

Message:

added function decimate_dem

Location:

inundation/ga/storm_surge/pyvolution

Files:

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

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -2385 +2385 @@
 
 
+def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
+                 verbose=False):
+    """Read Digitial Elevation model from the following NetCDF format (.dem)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
+    """
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, zeros, sum, reshape
+
+    root = basename_in
+    inname = root + '.dem'
+
+    #Open existing netcdf file to read
+    infile = NetCDFFile(inname, 'r')
+    if verbose: print 'Reading DEM from %s' %inname
+
+    #Read metadata
+    ncols = infile.ncols[0]
+    nrows = infile.nrows[0]
+    xllcorner = infile.xllcorner[0]
+    yllcorner = infile.yllcorner[0]
+    cellsize = infile.cellsize[0]
+    NODATA_value = infile.NODATA_value[0]
+    zone = infile.zone[0]
+    false_easting = infile.false_easting[0]
+    false_northing = infile.false_northing[0]
+    projection = infile.projection
+    datum = infile.datum
+    units = infile.units
+
+    dem_elevation = infile.variables['elevation']
+
+    assert len(dem_elevation) == nrows * ncols, 'elevation data do not match nrows*ncols'
+
+    #Get output file name
+    if basename_out == None:
+        outname = root + '_' + repr(cellsize_new) + '.dem'
+    else:
+        outname = basename_out + '.dem'
+
+    if verbose: print 'Write decimated NetCDF file to %s' %outname
+
+    #Determine some dimensions for decimated grid
+    (nrows_stencil, ncols_stencil) = stencil.shape
+    x_offset = ncols_stencil / 2
+    y_offset = nrows_stencil / 2
+    cellsize_ratio = cellsize_new / cellsize
+    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
+    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
+
+    #Open netcdf file for output
+    outfile = NetCDFFile(outname, 'w')
+
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF DEM format for compact and portable storage ' +\
+                           'of spatial point data'
+    #Georeferencing
+    outfile.zone = zone
+    outfile.projection = projection
+    outfile.datum = datum
+    outfile.units = units
+
+    outfile.cellsize = cellsize_new
+    outfile.NODATA_value = NODATA_value
+    outfile.false_easting = false_easting
+    outfile.false_northing = false_northing
+
+    outfile.xllcorner = xllcorner + (x_offset * cellsize)
+    outfile.yllcorner = yllcorner + (y_offset * cellsize)
+    outfile.ncols = ncols_new
+    outfile.nrows = nrows_new
+
+
+    # dimension definition
+    outfile.createDimension('number_of_points', nrows_new*ncols_new)
+
+    # variable definition
+    outfile.createVariable('elevation', Float, ('number_of_points',))
+
+    # Get handle to the variable
+    elevation = outfile.variables['elevation']
+
+    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
+
+    #Store data
+    global_index = 0
+    for i in range(nrows_new):
+        if verbose: print 'Processing row %d of %d' %(i, nrows)
+        lower_index = global_index
+        telev =  zeros(ncols_new, Float)
+        local_index = 0
+        trow = i * cellsize_ratio
+
+        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))
+            global_index += 1
+            local_index += 1
+
+        upper_index = global_index
+
+        elevation[lower_index:upper_index] = telev
+
+    assert global_index == nrows_new*ncols_new, 'index not equal to number of points'
+
+    infile.close()
+    outfile.close()
+

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -1819 +1819 @@
 
 
+    def test_decimate_dem(self):
+        """Test decimation of dem file
+        """
+
+        import os
+        from Numeric import ones, allclose, Float, arange
+        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
+
+        fid.ncols = ncols
+        fid.nrows = nrows
+        fid.xllcorner = 2000.5
+        fid.yllcorner = 3000.5
+        fid.cellsize = 25
+        fid.NODATA_value = -9999
+
+        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']
+
+        elevation[:] = (arange(nrows*ncols))
+
+        fid.close()
+
+        #generate the elevation values expected in the decimated file
+        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
+                         (  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,
+                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
+                         ( 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,
+                         (152+153+154+170+171+172+188+189+190) / 9.0,
+                         (156+157+158+174+175+176+192+193+194) / 9.0,
+                         (216+217+218+234+235+236+252+253+254) / 9.0,
+                         (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')
+
 
 #-------------------------------------------------------------
@@ -1824 +1907 @@
     #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')
     runner = unittest.TextTestRunner()
     runner.run(suite)