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

45

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):

Raster en niveaux de gris représentant l'offre de pollinisateurs sur un paysage

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

entrez la description de l'image ici

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.

James
la source
2
Qu'en est-il de ReadAsArray () et de la lecture dans le point? Alors seulement lire la cellule unique qui vous intéresse? Vous devez convertir les points coords en espace pixel et extraire la cellule nécessaire.
Jay Laura
1
Regardez le code de gdalsrsinfo, il vous montre comment utiliser GDALInvertGeoTransform () et basculer entre espace géographique et espace de pixels: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp
Si cela ne vous dérange pas d'utiliser PostGIS, voyez ceci . Il est extrêmement rapide et ne représente qu'une ligne SQL.
Mlt
Je garderai cela à l'esprit si je rencontre ce problème et si j'ai accès à une base de données PostGIS! Je ne l'ai pas fait pour ce problème particulier, alors la solution GDAL ci-dessous a fait l'affaire. Merci quand même!
James
@kyle Je ne sais pas si les choses ont changé, mais il semblerait que GDALInvGeoTransform ne soit pas inversé et il s'agit d'un exemple .
Mlt

Réponses:

61

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
utilisateur2856
la source
Et comment pouvez-vous enregistrer en tant que csv, table ou autre objet? Un objet avec la même longueur de coordonnées objetct
gonzalez.ivan90
Voici la question, Luke. Merci! gis.stackexchange.com/questions/269603/…
gonzalez.ivan90
les lignes qui assignent px/ pyont tort dans le cas où mx / my est en dehors des limites de rb, car int(-0.5) == 0. Vous avez besoin floor(...), puis vous devez vérifier qu'aucun des px/ pyn'est inférieur à zéro (ou le faire avant d'appeler int()), car les indices négatifs fonctionnent (ils obtiennent l'autre côté du tableau). J'aimerais savoir s'il existe un moyen plus ordonné de traiter ce problème. Comment réécrire ces lignes pour qu'elles traitent correctement les rotations?
naught101