[gdal-dev] basic raster math

[Bobby Braswell]

Hi Greg, here's a basic example without GDAL (and thus with no  
projection info). I hope it helps (and I hope there aren't typos).  --  
Rob Braswell


# load the modules
import sys
import numpy
import Image
from pylab import *

# open the files
im1 = Image.open("test1.tif")
im2 = Image.open("test2.tif")

# 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)

datadiff = data2 - data1

# set mode to interactive
matplotlib.interactive(True)

# plot a cross section
figure(1)
plot(datadiff[:,1000])

# look at the image
figure(2)
imshow(datadiff)



On Feb 14, 2008, at 4:24 PM, Greg Fiske wrote:

> Greetings,
>
> 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?
>
> Thanks,
>
> Greg
>
> Gregory Fiske
> Research Associate
> GIS and Remote Sensing Laboratory
> The Woods Hole Research Center
> 149 Woods Hole Road
> Falmouth, Massachusetts 02540
> 508-540-9900
> http://www.whrc.org
> gfiske at whrc.org
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

--
Bobby H. Braswell (Rob)
Institute for the Study of Earth Oceans and Space
University of New Hampshire
Durham, NH 03824 USA
1-603-862-2264