¿Obtener el valor de píxel del ráster GDAL bajo el punto OGR sin NumPy?

Estoy trabajando en un modelo computacional de la abundancia de polinizadores salvajes en un paisaje. El modelo en sí está completo, y ahora estoy luchando con un paso de procesamiento posterior.

Tengo mi ráster de suministro de polinizador GDAL que se parece a esto (los colores más claros significan una mayor visita de polinizador a un píxel):

Y tengo un archivo de forma OGR de puntos que representan ubicaciones de muestra en el paisaje:

Estoy tratando de ejecutar un análisis en los píxeles debajo de estos puntos, pero para hacerlo, necesito poder extraer el valor de un píxel debajo de un punto.

¿Es posible extraer el valor de un píxel debajo de un punto usando solo OGR y GDAL a través de Python? Preferiría evitar leer todo el ráster en la memoria ReadAsArray(), ya que mis rásteres de salida son muy, muy grandes (demasiado grandes para caber en la memoria).

Noté esta publicación , que es similar, pero requiere una llamada de línea de comandos.

Puede usar el método gdal.Dataset o gdal.Band ReadRaster. Consulte los tutoriales de la API de GDAL y OGR y el ejemplo a continuación. ReadRaster no utiliza / requiere numpy, el valor de retorno son datos binarios sin procesar y deben descomprimirse utilizando el módulo de estructura python estándar .

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

Alternativamente, dado que la razón que dio para no usar numpyfue evitar leer toda la matriz al usar ReadAsArray(), a continuación se muestra un ejemplo que usa numpyy no lee todo el ráster.

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