[gdal-dev] Getting data values from a first GeoTIFF and do some map algebra type operations with a second GeoTIFF

[jdmorgan]

Hello,

I am new to GDAL and raster programming in general and hoping that I can 
get some guidance. I have GeoTIFF which are projected in WGS84 and a 
single band.I am trying to start simply with the concept of getting out 
the data values from the GeoTIFF from one file so that I compare them to 
another.However, I am a bit stuck at the point of getting out data 
values from a first GeoTIFF and do some map algebra type operations with 
a second GeoTIFF.Here is what I have so far:

from osgeo import gdal, gdalconst

import os, sys, time

import numpy

ds = gdal.Open('C:\\Data\\GIS\\python\\gdal\\CMIP3-tb1_2021_2050.tif', 
gdal.GA_ReadOnly)

cols = ds.RasterXSize

rows = ds.RasterYSize

bands = ds.RasterCount

print cols, rows, bands

gdal.AllRegister()

gt = ds.GetGeoTransform()

print gt

xOrigin = gt[0]

yOrigin = gt[3]

pixWidth = gt[1]

print "pixWidth " + str(pixWidth)

pixHeight = gt[5]

print "pixHeight " + str(pixHeight)

for r in range(rows):

data = ds.GetRasterBand(1).ReadAsArray(0, r, cols, 1)

print str(data)

Is this the correct way to be going about this?Any guidance or best 
practices/algorithms would be greatly appreciated.

Thanks,

Derek

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20111123/0ae28142/attachment-0001.html