¿Cómo puedo obtener el valor del píxel de una GDAL raster en virtud de un OGR punto sin NumPy?

Question

Estoy trabajando en un modelo computacional de la abundancia de los polinizadores silvestres a través de un paisaje. El modelo en sí es completa, y yo ahora estoy luchando con un post-procesamiento de paso.

Tengo mi GDAL polinizador de suministro de trama que se ve algo como esto (los colores más claros media superior polinizador de las visitas a un pixel):

Y tengo una OGR shapefile de puntos que representan lugares de muestra en el paisaje:

Estoy tratando de ejecutar un análisis sobre los píxeles en estos puntos, pero para ello, necesito ser capaz de extraer el valor de un pixel en un punto.

Es posible extraer el valor de un pixel en un punto utilizando sólo OGR y GDAL a través de Python? Yo preferiría evitar la lectura de la totalidad de la trama en la memoria a través de la ReadAsArray(), como mi salida de rásteres son muy, muy grande (demasiado grande para caber en la memoria).

Me di cuenta de este post, que es similar, pero requiere de una llamada de línea de comandos.

Answer 1

Usted puede utilizar el gdal.Conjunto de datos o de gdal.Banda ReadRaster método. Ver el GDAL y OGR tutoriales de la API y el ejemplo de abajo. ReadRaster no uso/requieren numpy, el valor de retorno es binario de datos y necesita ser descomprimido con el estándar de python struct módulo.

Véase también el QGIS "Punto de Muestreo de la Herramienta de" plug-in para un GUI manera de hacer esto.

Un ejemplo:

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
    intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

    print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value

Como alternativa, desde la razón que dio para no usar numpy fue para evitar la lectura de toda la matriz en el uso de ReadAsArray(), a continuación se muestra un ejemplo que utiliza numpy y no leer toda la trama.

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename) 
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    intval=rb.ReadAsArray(px,py,1,1)
    print intval[0] #intval is a numpy array, length=1 as we only asked for 1 pixel value