Obtener todos los valores de un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS

En otros posts se ha usado el método ‘identify’ de QgsRasterDataProvider para acceder a valores individuales o a todos ellos de un ráster. Tiene la particularidad de que, al final, permite generar una capa de puntos con patrón de retícula si el muestreo tiene como objetivo todos los valores. Sin embargo, existe un método de QgsRasterDataProvider más eficiente: ‘block’. Este genera un objeto de la clase QgsRasterBlock que, a través de su método ‘value’, hace posible la obtención de todos ellos en un sólo paso o modificarlos mediante el método ‘setValue’.

El código siguiente se utiliza para seleccionar todos los valores de un ráster 3 x 3 y cambiar sólo aquellos que correpondan al valor 4 por el de 25 (imagen a continuación). Posteriormente, usando el módulo GDAL de Python, se imprimen estos cambios como un nuevo ráster.

from osgeo import gdal, osr
import numpy as np

layer = iface.activeLayer()

provider = layer.dataProvider()

extent = provider.extent()

rows = layer.height()
cols = layer.width()

xmin = extent.xMinimum()
ymax = extent.yMaximum()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

print rows, cols

block = provider.block(1, extent, cols, rows)

values = [ [] for i in range(rows) ]

for i in range(rows):
    for j in range(cols):
        if block.value(i,j) == 4:
            block.setValue(i,j,25)
            values[i].append(block.value(i,j))
        else:
            values[i].append(block.value(i,j))

raster = np.array(values)

geotransform = [xmin, xsize, 0, ymax, 0, -ysize]

# Create gtif file with rows and columns from parent raster 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/aleatorio_block.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Int16)

##writting output raster
dst_ds.GetRasterBand(1).WriteArray( raster )
 
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
 
# setting spatial reference of output raster 
epsg = layer.crs().postgisSrid()
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dst_ds.SetProjection( srs.ExportToWkt() )
 
#Close output raster dataset 
dst_ds = None

Después de ejecutado el código anterior en la Python Console de QGIS se obtiene:

aleatorio_block1

donde puede observarse, mediante el uso de Value Tool plugin, que los cambios han sido producidos tal como se esperaba.

Esta entrada fue publicada en PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

4 respuestas a Obtener todos los valores de un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS

  1. Pingback: Cortando (clipping) un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS | El Blog de José Guerrero

  2. Hola José:

    Estoy haciendo algunas pruebas desde un pluggin con la clase QgsRasterBlock, he cargado un raster de elevación con una sola banda de 1233 píxeles de ancho por 631 de alto. Cuando intento recorrer el valor de cada pixel QGIS se cierra con un “minivolcado”. He simplificado el código al máximo para comprobar que no estoy cometiendo errores, pero aún así me sigue sucediendo…

    block = dataProvider.block(1, selectedLayer.extent(), selectedLayer.width(), selectedLayer.height())
            
    for y2 in range(block.height()):                    
        for x2 in range(block.width()):
            value = block.value(x2, y2)
            message = "x: {}".format(x2) + " y: {}".format(y2) + " value: {}".format(value)
            logging.debug(message)
    

    Cuando miro el fichero de log veo lo siguiente:

    x: 631 y: 0 value: 0.0
    

    Sólo ha llegado hasta la columna 631 (no ha llegado a pasar a la segunda fila) y el valor que recupera es 0, cuando el NODATAVALUE es -999 y en las estadísticas el valor mínimo es de más de 1000… Parece como si block no llegara a estar cargado por completo, pero sin embargo block.isValid() es TRUE

    Se te ocurre que puede estar pasando?

    Muchas gracias

    • Este código funciona como se espera:

      selectedLayer = iface.activeLayer()
      dataProvider = selectedLayer.dataProvider()
      
      rows = selectedLayer.height()
      cols = selectedLayer.width()
      
      block = dataProvider.block(1, selectedLayer.extent(), cols, rows)
      
      for x in range(rows):
          for y in range(cols):
              value = block.value(x, y)
              message = "x: {}".format(x) + " y: {}".format(y) + " value: {}".format(value)
              print message
      

      Saludos

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s