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

[Joaquim Luis]

Derek,

These are the type of operations that I recommend GMT (built with GDAL 
off course). With GMT5 ii as simple as

grd2xyz yourGeotiff > outfile

this saves  x,y,z triplets. Add -ZTLa to write only the z values (many 
other options available)

Joaquim

> 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
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

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