[gdal-dev] basic raster math

[Christopher Barker]

Bobby Braswell wrote:
> Hi Greg, here's a basic example without GDAL (and thus with no 
> projection info). 

if you want to write it back out as a geotiff, GDAL will help there. 
Take a look at some of the utilities that come with GDAL that are 
written in Python -- essentially what you need to do is open the file, 
then get the rasterband with:

numpy_array = band.ReadAsArray()

then do math with the numpy_arrays, then write it back out with 
something like:

gdal.Band.WriteAsArray(numpy_array)


A note about the PIL code provided:
> # extract the data
> im1data = im.getdata()
> im2data = im.getdata()
> 
> # get the dimensions
> dim = im1.size
> dim2 = im2.size
> 
> # check dimensions
> if (dim != dim2):
>     print "images are not the same size"
>     sys.exit(0)
> 
> # convert the data from image type to array type
> # (there are a lot of implications to this but for now it should do what 
> you want)
> data1 = numpy.array(im1data).reshape(dim)
> data2 = numpy.array(im2data).reshape(dim)

All the above code can now be replaced with:

data1 = numpy.asarray(im1) (im1 is a PIL image object)

thanks to the "array interface"

-Chris



-- 
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov