Получить значение пикселя растра GDAL под точкой OGR без NumPy?

45

Я работаю над вычислительной моделью обилия диких опылителей в ландшафте. Сама модель завершена, и сейчас я борюсь с этапом постобработки.

У меня есть мой растровый растр GDAL, который выглядит примерно так (более светлые цвета означают большее посещение опылителя на пиксель):

Greyscale растр, представляющий поставку опылителей на ландшафте

И у меня есть шейп-файл OGR точек, представляющих образцы мест на ландшафте:

введите описание изображения здесь

Я пытаюсь провести некоторый анализ пикселей под этими точками, но для этого мне нужно иметь возможность извлечь значение пикселя под точкой.

Можно ли извлечь значение пикселя под точкой, используя только OGR и GDAL через Python? Я бы предпочел не читать весь растр в память ReadAsArray(), так как мои выходные растры очень, очень большие (слишком большие, чтобы поместиться в память).

Я заметил этот пост , который похож, но требует вызова из командной строки.

Джеймс
источник
2
Как насчет ReadAsArray () и только чтение в точку? Так что читайте только одну ячейку, которая вас интересует? Вам нужно будет преобразовать координаты точки в пиксельное пространство и извлечь необходимую ячейку.
Джей Лора
1
Посмотрите на код для gdalsrsinfo, он показывает, как использовать GDALInvertGeoTransform () и переключаться между географическим пространством и пространством пикселей: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp
Если вы не против использования PostGIS, посмотрите это . Это очень быстро и просто 1 строка SQL.
млт
Я буду помнить об этом, если столкнусь с этой проблемой и получу доступ к базе данных PostGIS! Я не для этой конкретной проблемы, поэтому решение GDAL ниже сделал свое дело. Спасибо хоть!
Джеймс
@kyle Я не знаю, если что-то изменилось, но похоже, что это GDALInvGeoTransform не инвертировать, и это пример .
млт

Ответы:

61

Вы можете использовать метод gdal.Dataset или gdal.Band ReadRaster. См. Учебники по GDAL и OGR API и пример ниже. ReadRaster не использует / требует numpy, возвращаемое значение является необработанными двоичными данными и должно быть распаковано с использованием стандартного модуля структуры python .

Пример:

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

В качестве альтернативы, поскольку причина, по которой вы отказались от использования, numpyзаключалась в том, чтобы не читать весь используемый массив ReadAsArray(), ниже приведен пример, в котором используется numpyи не читается весь растр.

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
user2856
источник
И как вы можете сохранить как CSV, таблицу или другой объект? Объект с той же длиной координат объекта
gonzalez.ivan90
Вот вопрос, Люк. Благодарность! gis.stackexchange.com/questions/269603/…
gonzalez.ivan90
строки, которые назначают px/ pyявляются неправильными в случае, когда mx / my лежит за пределами rb, потому что int(-0.5) == 0. Вам нужно floor(...), а затем вам нужно проверить, что ни один из px/ pyне меньше нуля (или сделать это перед вызовом int()), потому что работают отрицательные индексы (они получают другую сторону массива). Я хотел бы знать, есть ли более аккуратный способ справиться с этой проблемой. Кроме того, как вы переписываете эти строки, чтобы они правильно обрабатывали повороты?
naught101