Función en PyQGIS para crear un ráster de n filas x m columnas con values aleatorios

Para probar procedimientos y algoritmos sobre archivos ráster es conveniente que contengan values sencillos que nos permitan un mayor control sobre los cálculos y extensiones más pequeñas con el fin de reducir los tiempos de cómputo. Por ello, se piensa en la creación de rásteres síntéticos.

Para producir un ráster sintético lo más adecuado es que esté alineado con el ráster base y ello se logra copiando los parámetros importantes del original (vértice de inicio, resolución y proyección) en el ráster destino. Para ello, nos haremos servir de las funciones propias de GDAL; aunque el provider de PyQGIS es el que nos permitirá determinar el directorio necesario para poder acceder a la dataset del ráster que corresponde a la capa activa.

Por otra parte, los values aleatorios del ráster son generados a través de numpy, con la dimensión que dicte nuestro deseo, para ser escritos al destinatario utilizando otra vez GDAL. El código, en forma de función, se resume a continuación:

#-*- coding: utf-8 -*-

from qgis.gui import QgsMessageBar
import numpy as np
from osgeo import osr, gdal
from gdalconst import *
from os import path

def rasterCloneLayerRand(rows, columns, v_min = 1, v_max = 10):

    layer = iface.activeLayer()
    
    if layer is None or layer.type() != QgsMapLayer.RasterLayer:
        
        iface.messageBar().pushMessage("Warning:",
                                       "There is not a Raster Layer",
                                       QgsMessageBar.WARNING, 5)
        
        return

    provider = layer.dataProvider()
     
    my_path= provider.dataSourceUri()
     
    (root, filename) = path.split(my_path)
     
    dataset = gdal.Open(my_path)

    prj = dataset.GetProjection()
     
    band = dataset.GetRasterBand(1)

    #Se obtiene la metadata del ráster
     
    #geotransform[0] = top left x
    #geotransform[1] = w-e pixel resolution
    #geotransform[3] = top left y
    #geotransform[5] = n-s pixel resolution
     
    geotransform = dataset.GetGeoTransform()

    # Establece el nombre del ráster de salida
    output_file = root + "/" + "out2.tif"

    # Crea el gtif con las filas y columnas del ráster padre 
    driver = gdal.GetDriverByName("GTiff")
    
    dst_ds = driver.Create(output_file, 
                           columns, 
                           rows, 
                           1, 
                           band.DataType)

    #Declara el arreglo del raster 
    raster = np.random.randint(v_min, 
                               v_max + 1, 
                               size = (rows, columns))

    ##Escribe al ráster de salida
    dst_ds.GetRasterBand(1).WriteArray( raster )

    #Establece la extensión del ráster de salida
    # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
    dst_ds.SetGeoTransform(geotransform)
      
    # Establece la referencia espacial del ráster de salida 
    srs = osr.SpatialReference(wkt = prj)
    dst_ds.SetProjection( srs.ExportToWkt() )

    ##Cierra el dataset del ráster padre
    dataset = None

    #Cierra el dataset del ráster de salida
    dst_ds = None

Cuando ejecutamos la función por defecto en la Python Console, se le pasan como parámetros el número de filas y el número de columnas para producir rásteres con valores enteros entre 1 y 10; ambos inclusive. Si se quieren valores diferentes entonces hay que suministrar adicionalmente el valor mínimo y el máximo, es decir, en total cuatro parámetros.

Para obtener un ráster de 30 x 30 alineado con el ráster base utah_demUTM2.tif (ver figura), con values entre 1 y 3, se tiene que rasterCloneLayerRand(30,30,1,3).

random1

Para producir un ráster de 30 x 30, con los valores por defecto entre 1 y 10, tenemos que rasterCloneLayerRand(30,60).

random2

El plugin Value Tool en la imagen anterior nos permite corroborar que los ráster fueron producidos de la manera prevista.

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

Una respuesta a Función en PyQGIS para crear un ráster de n filas x m columnas con values aleatorios

  1. Pingback: Cómo introducir valores numéricos en un plugin de QGIS a través de una QLineEdit | El Blog de José Guerrero

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