[Gdal-dev] Numpy and GDAL

[Frank Warmerdam]

Folks,

I haven't really mastered the zen of Numpy (aka Numeric), the array math
package for Python, but I did a nice hack today that I thought might be of
interest to some others.  The requirement was to add "nodata" support to
the gdal_merge.py script, a simple script to merge a bunch of raster files
into one mosaic.

This devolves down to a function call to copy a rectangle of data from a
source file to a destination file.  I changed it to look like this:

# =============================================================================
def raster_copy_with_nodata( s_fh, s_xoff, s_yoff, s_xsize, s_ysize, s_band_n,
                              t_fh, t_xoff, t_yoff, t_xsize, t_ysize, t_band_n,
                              nodata ):

     if verbose != 0:
         print 'Copy %d,%d,%d,%d to %d,%d,%d,%d.' \
               % (s_xoff, s_yoff, s_xsize, s_ysize,
              t_xoff, t_yoff, t_xsize, t_ysize )

     s_band = s_fh.GetRasterBand( s_band_n )
     t_band = t_fh.GetRasterBand( t_band_n )

     data_src = s_band.ReadAsArray( s_xoff, s_yoff, s_xsize, s_ysize,
                                    t_xsize, t_ysize )
     data_dst = t_band.ReadAsArray( t_xoff, t_yoff, t_xsize, t_ysize )

     nodata_test = Numeric.equal(data_src,nodata)
     to_write = Numeric.choose( nodata_test, (data_src, data_dst) )

     t_band.WriteArray( to_write, t_xoff, t_yoff )

     return 0

The key items are:
  o Read source window as an "array" at the same resolution as the output
    window.
  o Read the existing data in the destination file for the target window
    as an array.
  o Create an array of 0/1 values indicating which pixels in the source
    array match the nodata value - "Numeric.equal(data_src,nodata)" where
    nodata is a simple scalar.
  o Use the Numeric.choose() function to choose the existing target file
    value where the test returned 1, and the source data where it returned 0.
  o Write to destination file.

The nice thing is that Numeric.equal() and Numeric.choose() are internally
implemented in C so there is no per-pixel python script evaluation going on.

Note that the Numeric.choose() function can take a longer vector of values.
For instance, you could have a 256 entry lut and use an 8bit input image
array to do a lookup in the lut.

Pretty nice stuff!

PS. I have committed a change to gdal_merge.py so that it now supports
"nodata" values ... pixel values in the source image that should not be
copied to the merged image.

Later,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent