Obtenir la valeur en pixels du raster GDAL sous le point OGR sans NumPy?

Je travaille sur un modèle informatique de l'abondance de pollinisateurs sauvages à travers un paysage. Le modèle lui-même est complet et je suis maintenant aux prises avec une étape de post-traitement.

J'ai mon raster d'approvisionnement en pollinisateur GDAL qui ressemble à ceci (des couleurs plus claires signifient une fréquentation plus élevée des pollinisateurs au pixel près):

Et j'ai un fichier de formes OGR de points représentant des exemples d'emplacements dans le paysage:

J'essaie d'analyser les pixels situés sous ces points, mais pour ce faire, je dois pouvoir extraire la valeur d'un pixel sous un point.

Est-il possible d'extraire la valeur d'un pixel sous un point en utilisant uniquement OGR et GDAL via Python? Je préférerais éviter de lire l'intégralité du raster en mémoire ReadAsArray()car mes rasters en sortie sont très très gros (trop gros pour être stockés dans la mémoire).

J'ai remarqué ce message , qui est similaire, mais nécessite un appel de ligne de commande.

Vous pouvez utiliser la méthode ReadRaster gdal.Dataset ou gdal.Band . Voir les tutoriels sur les API GDAL et OGR et l'exemple ci-dessous. ReadRaster n'utilise pas / require numpy, la valeur de retour correspond à des données binaires brutes et doit être décompressée à l'aide du module struct Python standard .

Un exemple:

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

Alternativement, puisque la raison que vous avez indiquée pour ne pas utiliser numpyétait d'éviter de lire le tableau entier en utilisant ReadAsArray(), voici un exemple qui utilise numpyet ne lit pas le raster en entier.

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