Uzyskujesz wartość piksela rastra GDAL poniżej punktu OGR bez NumPy?

Pracuję nad obliczeniowym modelem liczebności dzikich zapylaczy w krajobrazie. Sam model jest kompletny, a teraz mam problem z etapem przetwarzania końcowego.

Mam raster zapylacza GDAL, który wygląda mniej więcej tak (jaśniejsze kolory oznaczają większą wizytę zapylacza na piksel):

I mam plik kształtu OGR punktów reprezentujących przykładowe lokalizacje w krajobrazie:

Próbuję przeprowadzić analizę pikseli pod tymi punktami, ale aby to zrobić, muszę być w stanie wyodrębnić wartość piksela poniżej punktu.

Czy jest możliwe wyodrębnienie wartości piksela pod punktem przy użyciu tylko OGR i GDAL przez Python? Wolałbym unikać wczytywania całego rastra do pamięci ReadAsArray(), ponieważ moje wyjściowe rastry są bardzo, bardzo duże (zbyt duże, aby zmieściły się w pamięci).

Zauważyłem ten post , który jest podobny, ale wymaga wywołania z wiersza poleceń.

Możesz użyć metody gdal.Dataset lub gdal.Band ReadRaster. Zobacz samouczki interfejsu API GDAL i OGR oraz poniższy przykład. ReadRaster nie używa / wymaga numpy, zwracaną wartością są surowe dane binarne i należy je rozpakować przy użyciu standardowego modułu strukturalnego Pythona .

Przykład:

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

Alternatywnie, ponieważ powodem, dla którego podałeś, że nie używasz, numpybyło uniknięcie odczytu całej tablicy podczas używania ReadAsArray(), poniżej znajduje się przykład, który używa numpyi nie czyta całego rastra.

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