Pixelwert des GDAL-Rasters unter OGR-Punkt ohne NumPy erhalten?

Ich arbeite an einem Rechenmodell für die Häufigkeit von Wildbestäubern in einer Landschaft. Das Modell selbst ist vollständig, und ich habe jetzt mit einem Nachbearbeitungsschritt zu kämpfen.

Ich habe mein GDAL-Bestäuberversorgungs-Raster, das ungefähr so ​​aussieht (hellere Farben bedeuten einen höheren Bestäuberaufruf pro Pixel):

Und ich habe ein OGR-Shapefile mit Punkten, die Beispielpositionen in der Landschaft darstellen:

Ich versuche, einige Analysen an den Pixeln unter diesen Punkten durchzuführen, aber dazu muss ich in der Lage sein, den Wert eines Pixels unter einem Punkt zu extrahieren.

Ist es möglich, den Wert eines Pixels unter einem Punkt nur mit OGR und GDAL über Python zu extrahieren? Ich würde es vorziehen, das Lesen des gesamten Rasters in den Speicher zu vermeiden ReadAsArray(), da meine Ausgabe-Raster sehr, sehr groß sind (zu groß, um in den Speicher zu passen).

Mir ist dieser Beitrag aufgefallen , der ähnlich ist, aber einen Befehlszeilenaufruf erfordert.

Sie können die Methode gdal.Dataset oder gdal.Band ReadRaster verwenden. Weitere Informationen finden Sie in den Tutorials zur GDAL- und OGR- API sowie im folgenden Beispiel. ReadRaster verwendet / erfordert kein numpy, der Rückgabewert besteht aus rohen Binärdaten und muss mit dem Standard-Python- Strukturmodul entpackt werden .

Ein Beispiel:

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

Da Sie als Grund für die Nichtverwendung angegeben haben, dass nicht numpydas gesamte Array in using gelesen werden soll ReadAsArray(), finden Sie im Folgenden ein Beispiel, das numpydas gesamte Raster verwendet und nicht einliest.

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