Cómo modificar individualmente los valores de un ráster en QGIS con Python

Para modificar individualmente los valores de un ráster hay que acceder al array. En PyQGIS es factible la modificación mediante el método ‘write’ de QgsRasterDataProvider; tal como señalan aquí. Sin embargo, no he podido encontrar un ejemplo práctico de cómo se usa el puntero a través de la clase QgsRasterPipe. Por tanto, he apelado a los módulos gdal y struct de Python, a través de una línea de barrido (scanline), para obtener la matriz de valores del ráster (lista de tuplas).

Como las tuplas son inmutables, hay que hacer una conversión previa a listas para modificar los valores individuales por (fila, columna). Una vez hecha la sustitución se convierte la lista de listas en array para que pueda ser incorporada en el nuevo ráster a crear preservando así la integridad del original.

El código siguiente realiza la sustitución en la (fila, columna) = (1, 5) por el valor 15.

from osgeo import gdal, osr
import os, struct
import numpy as np

def changeRasterValue(band, row, column, value):
    
    fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
    
    data_type = band.DataType
    
    BandType = gdal.GetDataTypeName(band.DataType)

    raster = []
        
    for y in range(band.YSize):

        scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, data_type)
        values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
        raster.append(values)

    raster = [ list(item) for item in raster ]
    
    raster[row][column] = value

    #transforming list in array
    raster = np.asarray(np.reshape(raster, (band.YSize, band.XSize)))
    
    return raster

layer = iface.activeLayer()

provider = layer.dataProvider()

path = provider.dataSourceUri()

(raiz, filename) = os.path.split(path)

dataset = gdal.Open(path)

#Get projection
prj = dataset.GetProjection()

#setting band
number_band = 1

band = dataset.GetRasterBand(number_band)

#Get raster metadata
geotransform = dataset.GetGeoTransform()

# Set name of output raster
output_file = "/home/zeito/pyqgis_data/neighbors_output2.tif"
 
# Create gtif file with rows and columns from parent raster 
driver = gdal.GetDriverByName("GTiff")

row = 1 
column = 5
value = 15

raster = changeRasterValue(band, row, column, value)
    
dst_ds = driver.Create(output_file, 
                       band.XSize, 
                       band.YSize, 
                       number_band, 
                       band.DataType)

#writting output raster
dst_ds.GetRasterBand(number_band).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 
srs = osr.SpatialReference(wkt = prj)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

#Close main raster dataset
dataset = None

Usando un ráster aleatorio (21 filas x 21 columnas) con valores de 1 a 10, después de correr el script, la sustitución fue satisfactoria en la posición seleccionada del ráster de salida; tal como se puede observar en la imagen siguiente:

changevalues1

usando el plugin ‘Value Tool’ para corroborarlo.

Esta entrada fue publicada en Código Python, PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

4 respuestas a Cómo modificar individualmente los valores de un ráster en QGIS con Python

  1. Will Quenta dijo:

    hay alguna herramienta que pueda hacer la reclasificacion de una raster solo por el area de interes como lo hace el recode de ERDAS tengo una raster de gran extension y se confudieron algunos valores y solo quiero cambiar esos valores sin crear un nuevo archivo de salida

  2. Will Quenta dijo:

    Gracias, te lo agradezco, pero tengo el siguiente problema. Tengo un ráster que es producto de una reclasificación y al utilizar el ‘Value Tool’ sólo me muestra los values y las coordenadas. Cómo hago para que me muestre fila,columna?

  3. No es posible sino modificas por tu cuenta el plugin. Para ello hay que modificar dos formas *.py del mismo. Primero, tienes que pasarle una referencia de la posición y de las capas a la función ‘printInTable’ de ‘valuewidget.py’ y después transformar la posición en índices de fila, columna usando las coordenadas del vértice superior y la resolución de cada ráster. Segundo, tienes que modificar la forma que contiene los objetos Qt para que incorpore los valores de (fila, columna); además de la layer y el value.

    Se vería de la forma que está en la imagen siguiente:

    value_tool3

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