Mendapatkan nilai piksel raster GDAL di bawah titik OGR tanpa NumPy?

Saya sedang mengerjakan model komputasi dari banyaknya penyerbuk liar di seluruh lanskap. Model itu sendiri selesai, dan saya sekarang berjuang dengan langkah post-processing.

Saya memiliki raster pasokan penyerbuk GDAL saya yang terlihat seperti ini (warna yang lebih terang berarti kunjungan penyerbuk yang lebih tinggi ke satu piksel):

Dan saya memiliki OGR shapefile poin yang mewakili lokasi sampel di lanskap:

Saya mencoba menjalankan beberapa analisis pada piksel di bawah titik-titik ini, tetapi untuk melakukannya, saya harus dapat mengekstraksi nilai piksel di bawah satu titik.

Apakah mungkin untuk mengekstraksi nilai piksel di bawah satu titik menggunakan hanya OGR dan GDAL melalui Python? Saya lebih suka menghindari membaca seluruh raster ke dalam memori ReadAsArray(), karena raster output saya sangat, sangat besar (terlalu besar untuk masuk ke dalam memori).

Saya perhatikan posting ini , yang mirip, tetapi membutuhkan panggilan baris perintah.

Anda dapat menggunakan metode ReadRaster gdal.Dataset atau gdal.Band . Lihat tutorial GDAL dan OGR API dan contoh di bawah ini. ReadRaster tidak menggunakan / membutuhkan numpy, nilai kembali adalah data biner mentah dan perlu dibongkar menggunakan modul struct python standar .

Sebuah contoh:

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

Atau, karena alasan yang Anda berikan untuk tidak menggunakan numpyadalah untuk menghindari membaca seluruh array yang digunakan ReadAsArray(), di bawah ini adalah contoh yang menggunakan numpydan tidak membaca seluruh raster di.

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