Raster Query com a GDAL

Problema: saber o valor do pixel de uma imagem georeferenciada para um determinado par de coordenadas e uma determinada banda.

Solução utilizando GDAL Python bindings:

Nota: no código que se segue substituir underscore(_) por espaço em branco.

import struct
import gdal
from gdalconst import *

filename = "image.tif"
dataset = gdal.Open( filename, GA_ReadOnly )

def CellValue(dataset, band, x, y):
_band = dataset.GetRasterBand(band)
_if (band.GetNoDataValue() == None):
___band.SetNoDataValue(-9999)
_ndv = band.GetNoDataValue()
_cols = band.XSize
_rows = band.YSize
_geotransform = dataset.GetGeoTransform()
_cellSizeX = geotransform[1]
_cellSizeY = -1 * geotransform[5]
_minx = geotransform[0]
_maxy = geotransform[3]
_maxx = minx + (cols * cellSizeX)
_miny = maxy - (rows * cellSizeY)
_print 'bbox(real-world coords):',minx,',',miny,',',maxx,',',maxy
_if ((x <> maxx) or (y <> maxy)):
___print 'given point does not fall within grid'
___return ndv
_xLoc = (x - minx) / cellSizeX
_xLoc = int(xLoc)
_yLoc = (maxy - y) / cellSizeY
_yLoc = int(yLoc)
_#print 'point (pixels): ',xLoc,',',yLoc
_if ((xLoc <> cols - 0.5)):
___#print 'xcoordinate out of bounds'
___return ndv
_if ((yLoc <> rows - 0.5)):
___#print 'ycoordinate out of bounds'
___return ndv
_strRaster = band.ReadRaster(xLoc, yLoc, 1, 1, 1, 1, band.DataType)
_sDT = gdal.GetDataTypeName(band.DataType)
_if (sDT == 'Int16'):
___dblValue = struct.unpack('h', strRaster)
_elif (sDT == 'Float32'):
___dblValue = struct.unpack('f', strRaster)
_elif (sDT == 'Byte'):
___dblValue = struct.unpack('B', strRaster)
_else:
___print 'unrecognized DataType:', gdal.GetDataTypeName(band.DataType)
___return ndv
_print sDT
_return dblValue[0]

Testar a função:

print CellValue(dataset, 1, 222960.0, 285040.0)

Fonte: thread da mailing-list gdal-dev.