[gdal-dev] basic raster math

[Jose Luis Gomez Dans]

Hi!
> 
> I'm a new python user, so pardon my simple question.  I'd like to
> do some simple raster math within Python.  For example, I'd like
> to read in 2 Geotiffs and add one to the other -- creating a new
> output Geotiff.  I've experimented with the gdal and Numpy
> libraries, but I haven't made much progress on my own.  Can
> somebody send me some example code?

I do this all the time, and I am yet to find a better combination than Python + gdal, so that's good news!

To read data in (I am still using the old generation Python bindings):
import gdal
#fname is the file you want to read in, e.g. fname="thingy.tif"
g = gdal.Open(fname)
data = numpy.array(g.GetRasterBand(1).ReadAsArray().tolist())
#The previous line converts the oldgen python Numeric array into a 
#spanking numpy array.
data= numpy.where(numpy.isnan(data),0.0, data)#... which you can use
#normally here.

In order to write a file, I have a function that just does that. Note that I do analyses on overlapping regions, with the same projections and so on and so forth, so it is specialised, but you can tailor it to your needs:

def WriteRaster (dst_filename, raster):
	import osr
	format = "GTiff"
	driver = gdal.GetDriverByName( format )
	dst_ds = driver.Create( dst_filename, 71, 73,\
               1,gdal.GDT_Float32,options=["COMPRESS=PACKBITS","TFW=YES"] )
	dst_ds.SetGeoTransform( [-19.5, 1.0, 0.0, 37.5, 0.0, -1.0] )
	srs = osr.SpatialReference()
	srs.ImportFromEPSG(4326) #WGS84 lat long.
	dst_ds.SetProjection( srs.ExportToWkt() )
	dst_ds.GetRasterBand(1).WriteArray( raster )
	dst_ds = None

Hope that helps!
Jose
-- 
Der GMX SmartSurfer hilft bis zu 70% Ihrer Onlinekosten zu sparen! 
Ideal für Modem und ISDN: http://www.gmx.net/de/go/smartsurfer