Obtendo o valor de pixel da varredura GDAL sob o ponto OGR sem NumPy?

45

Estou trabalhando em um modelo computacional da abundância de polinizadores selvagens em uma paisagem. O modelo em si está completo e agora estou lutando com uma etapa de pós-processamento.

Eu tenho a minha varredura de suprimento de polinizadores GDAL que se parece com isso (cores mais claras significam maior visitação do polinizador em um pixel):

Raster em escala de cinza representando o fornecimento de polinizadores em uma paisagem

E eu tenho um shapefile de pontos do OGR representando locais de amostra na paisagem:

insira a descrição da imagem aqui

Estou tentando executar algumas análises nos pixels sob esses pontos, mas, para fazer isso, preciso extrair o valor de um pixel sob um ponto.

É possível extrair o valor de um pixel em um ponto usando apenas OGR e GDAL através do Python? Eu preferiria evitar ler toda a varredura na memória ReadAsArray(), pois meus rasters de saída são muito, muito grandes (muito grandes para caber na memória).

Notei este post , que é semelhante, mas requer uma chamada de linha de comando.

James
fonte
2
E o ReadAsArray () e apenas a leitura no ponto? Então, leia apenas a única célula em que você está interessado? Você precisaria converter as cordas de ponto em espaço de pixel e extrair a célula necessária.
Jay Laura
1
Veja o código para gdalsrsinfo, ele mostra como usar GDALInvertGeoTransform () e alternar entre espaço geográfico e pixel: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp
Se você não se importa em usar o PostGIS, veja isso . É extremamente rápido e é apenas uma linha SQL.
mlt
Vou ter isso em mente se me deparar com esse problema e tiver acesso a um banco de dados PostGIS! Como não resolvi esse problema em particular, a solução GDAL abaixo fez o truque. Obrigado, no entanto!
James
@kyle Eu não sei se as coisas mudaram, mas parece que ele é GDALInvGeoTransform não invertido e este é um exemplo .
Mlt 19/09

Respostas:

61

Você pode usar o método gdal.Dataset ou gdal.Band ReadRaster. Veja os tutoriais da API GDAL e OGR e o exemplo abaixo. O ReadRaster não usa / requer numpy, o valor de retorno são dados binários brutos e precisa ser descompactado usando o módulo de estrutura python padrão .

Um exemplo:

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

Como alternativa, como o motivo que você deu para não usar numpyfoi para evitar a leitura de toda a matriz em uso ReadAsArray(), abaixo está um exemplo que usa numpye não lê toda a varredura.

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
fonte
E como você pode salvar como CSV, tabela ou outro objeto? Um objeto com o mesmo comprimento de coordenadas objetct
gonzalez.ivan90
Aqui está a pergunta, Luke. Obrigado! gis.stackexchange.com/questions/269603/...
gonzalez.ivan90
as linhas que atribuem px/ pyestão erradas no caso em que mx / my está fora dos limites de rb, porque int(-0.5) == 0. Você precisa floor(...)e, em seguida, precisa verificar se nenhum de px/ pyé menor que zero (ou o faz antes de chamar int()), porque os índices negativos funcionam (eles ficam do outro lado da matriz). Gostaria muito de saber se existe uma maneira mais clara de lidar com esse problema. Além disso, como você reescreve essas linhas para que elas lidem corretamente com as rotações?
naught101 7/03