Cortando (clipping) un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS

En el post anterior se desarrolló un script para seleccionar todos los valores de un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS. Cualquier sub región de dicho ráster puede ser establecida simplemente con indicar su vértice superior izquierdo y su vértice inferior derecho. Por tanto, es lógico suponer que, una vez seleccionada dicha región, pueda ser escrita como un nuevo ráster. Esto sería equivalente a lo que produce la herramienta ‘Clipper’.

Para probarlo, se va a considerar la situación siguiente donde un ráster de 20 x 20 va ser “cortado” con base en los puntos que se observan en la imagen:

aleatorio1

La base del proceso consiste en averiguar los índices de fila y columna del ráster original que corresponden a cada punto y usarlos como base para encontrar los nuevos parámetros de ‘geotransform’ del raster “cortado”. Posteriormente, se clona el ‘provider’, a través de un objeto de la clase QgsRasterPipe, para escribir el ráster a través del método ‘writeRaster’ de QgsRasterFileWriter. El código completo se incluye a continuación:

#active layer
layer = iface.activeLayer()
#raster provider
provider = layer.dataProvider()
#original raster extent
extent = provider.extent()
#original raster resolution
rows = layer.height()
cols = layer.width()
#pixel size (original raster)
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()
#origin (original raster)
xmin = extent.xMinimum()
ymax = extent.yMaximum()

#cut points
pt1 = QgsPoint(355377.753377, 4473169.25631) #(3,5)
pt2 = QgsPoint(355682.174064, 4472867.27099) #(7,9)

points = [pt1, pt2]

#indexes for cut points
indxs = []

for point in points:
    new_row = int((ymax - point[1])/ysize)
    new_col = int((point[0] - xmin)/xsize)
    indxs.append([new_row, new_col])
    
new_pt1 = QgsPoint(xmin + indxs[0][1]*xsize, ymax - indxs[0][0]*ysize)
new_pt2 = QgsPoint(xmin + (indxs[1][1]+1)*xsize, ymax - (indxs[1][0]+1)*ysize)

new_rows = indxs[1][0] - indxs[0][0] + 1
new_cols = indxs[1][1] - indxs[0][1] + 1

pipe = QgsRasterPipe()
 
pipe.set(provider.clone())
 
rasterWriter = QgsRasterFileWriter("/home/zeito/pyqgis_data/aleatorio_write.tif")
 
CRS = layer.crs()

new_extent = QgsRectangle(new_pt1, new_pt2)
 
error = rasterWriter.writeRaster(pipe, new_cols, new_rows, new_extent, CRS)
 
if error == QgsRasterFileWriter.NoError:
    print "Raster was saved successfully!"
 
else:
    print "Raster was not saved!"

Después de ejecutado en la Python Console de QGIS se obtiene lo siguiente:

aleatorio2

El desfase ligero entre ambos rásteres (original y “cortado”), que se observa al hacer ‘Zoom In’, es menor que 1×10-3 metros por lo que el procedimiento se puede considerar exitoso.

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

Una respuesta a Cortando (clipping) un ráster mediante un objeto de la clase QgsRasterBlock en PyQGIS

  1. 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

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