¿Obtener el valor de píxel del ráster GDAL bajo el punto OGR sin NumPy?

45

Estoy trabajando en un modelo computacional de la abundancia de polinizadores salvajes en un paisaje. El modelo en sí está completo, y ahora estoy luchando con un paso de procesamiento posterior.

Tengo mi ráster de suministro de polinizador GDAL que se parece a esto (los colores más claros significan una mayor visita de polinizador a un píxel):

Ráster en escala de grises que representa el suministro de polinizadores en un paisaje

Y tengo un archivo de forma OGR de puntos que representan ubicaciones de muestra en el paisaje:

ingrese la descripción de la imagen aquí

Estoy tratando de ejecutar un análisis en los píxeles debajo de estos puntos, pero para hacerlo, necesito poder extraer el valor de un píxel debajo de un punto.

¿Es posible extraer el valor de un píxel debajo de un punto usando solo OGR y GDAL a través de Python? Preferiría evitar leer todo el ráster en la memoria ReadAsArray(), ya que mis rásteres de salida son muy, muy grandes (demasiado grandes para caber en la memoria).

Noté esta publicación , que es similar, pero requiere una llamada de línea de comandos.

James
fuente
2
¿Qué pasa con ReadAsArray () y solo leer en el punto? Entonces, ¿solo lee la celda única que le interesa? Necesitaría convertir de las coordenadas de punto a espacio de píxeles y extraer la celda necesaria.
Jay Laura
1
Mire el código de gdalsrsinfo, le muestra cómo usar GDALInvertGeoTransform () y cambie entre el espacio geográfico y el espacio de píxeles: trac.osgeo.org/gdal/browser/trunk/gdal/apps/gdalsrsinfo.cpp
Si no te importa usar PostGIS, mira esto . Es extremadamente rápido y es solo 1 línea SQL.
mlt
¡Lo tendré en cuenta si me encuentro con este problema y tengo acceso a una base de datos PostGIS! No lo hice para este problema en particular, por lo que la solución GDAL a continuación hizo el truco. Gracias, sin embargo!
James
@kyle No sé si las cosas han cambiado, pero parece que GDALInvGeoTransform no se invierte y este es un ejemplo .
mlt

Respuestas:

61

Puede usar el método gdal.Dataset o gdal.Band ReadRaster. Consulte los tutoriales de la API de GDAL y OGR y el ejemplo a continuación. ReadRaster no utiliza / requiere numpy, el valor de retorno son datos binarios sin procesar y deben descomprimirse utilizando el módulo de estructura python estándar .

Un ejemplo:

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

Alternativamente, dado que la razón que dio para no usar numpyfue evitar leer toda la matriz al usar ReadAsArray(), a continuación se muestra un ejemplo que usa numpyy no lee todo el ráster.

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
usuario2856
fuente
¿Y cómo puede guardar como csv, tabla u otro objeto? Un objeto con la misma longitud de coordenadas objetct
gonzalez.ivan90
Aquí está la pregunta, Luke. ¡Gracias! gis.stackexchange.com/questions/269603/…
gonzalez.ivan90
las líneas que asignan px/ pyestán equivocadas en el caso de que mx / my se encuentre fuera de los límites de rb, porque int(-0.5) == 0. Necesita floor(...), y luego debe verificar que ninguno de px/ pysea ​​menor que cero (o hacerlo antes de llamar int()), porque los índices negativos funcionan (obtienen el otro lado de la matriz). Me encantaría saber si hay una forma más ordenada de lidiar con este problema. Además, ¿cómo reescribe esas líneas para que se ocupen correctamente de las rotaciones?
naught101