Obtención del valor de píxel de ráster GDAL bajo el punto OGR sin NumPy?

39

Estoy trabajando en un modelo computacional de la abundancia de polinizadores silvestres 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 del polinizador a un píxel):

YtengounarchivodeformaOGRdepuntosquerepresentanubicacionesdemuestraenelpaisaje:

Estoy intentando realizar 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 en un punto utilizando solo OGR y GDAL a través de Python? Preferiría evitar leer todo el ráster en la memoria a través de 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 línea de comandos llamada.

    
pregunta James 11.01.2013 - 22:19

1 respuesta

56

Puede utilizar gdal.Dataset o gdal.Band método ReadRaster. Consulte GDAL y OGR tutoriales de API y el siguiente ejemplo. ReadRaster no usa / requiere números, el valor de retorno son datos binarios sin procesar y se debe desempaquetar usando el python estándar struct módulo.

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

Como alternativa, dado que la razón que dio para no usar numpy fue evitar la lectura de la matriz completa al usar ReadAsArray() , a continuación se muestra un ejemplo que usa numpy y no lee el ráster completo.

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
    
respondido por el Luke 12.01.2013 - 00:34

Lea otras preguntas en las etiquetas