Cómo recortar todos los ficheros ráster de una misma carpeta con PyQGIS

En el ejercicio de hoy vamos a ver un ejemplo de cómo recortar, utilizando la herramienta Clip ráster by mask, todos los ficheros ráster que tengamos almacenados en una carpeta concreta, automatizando el proceso con PyQGIS.

Para el desarrollo del ejercicio vamos a utilizar los ficheros ráster en formato *tif de la Red de Información Ambiental de Andalucía (REDIAM). En concreto, vamos a descargar los ficheros ráster de precipitación anual desde 1940 hasta 2019 para toda Andalucía. El enlace para descargar la carpeta completa.

Como ejemplo voy a seleccionar solamente los ficheros desde 2010, aunque si quisiéramos recortar todos los ráster descargados tan solo había que colocarlos en la misma carpeta. Para almacenar los ficheros de ejemplo he creado una carpeta llamada “raster” dentro de nuestro directorio de trabajo.

Para recortar los ráster vamos a utilizar una única máscara, que va a ser en este caso un fichero vectorial de la provincia de Cádiz. Para generarlo, podemos descargar del Centro de Descargas del IGN (CNIG), dentro de “Información Geográfica de Referencia”, los recintos provinciales. De esa capa extraeremos mediante una selección por atributos la provincia de Cádiz. Exportaremos la capa (Export – Save Features As) con el nombre “cadiz.shp” dentro de nuestro directorio de trabajo.

Para trabajar el código vamos a utilizar para mayor comodidad la consola de Python de QGIS Para acceder a ella: “Ctrl + Alt + P” y hacer click en el icono señalado.

Creamos un nuevo fichero .py y añadimos en primer lugar las librerías que nos van a hacer falta para ejecutar la herramienta, acceder a los directorios y crear carpetas.

/**
from osgeo import gdal
import os
from pathlib import Path
*/

Nombramos una nueva variable “basepath” que contendrá el directorio principal sobre el que vamos a trabajar. En mi caso, completo la ruta de la carpeta de usuario (“Path().home()”) que tengo definida añadiéndole con el método .joinpath() la ruta sobre la que quiero trabajar (“carpeta de usuario\\Documents\\entradas\\preci”).

Tras definir el directorio, creo las variables que van a contener los ficheros ráster (en este caso la carpeta “raster” que hemos creado dentro de la ruta del “basepath”) y el fichero vectorial (dentro del “basepath”, el fichero “cadiz.shp”).

/**
basepath = Path().home().joinpath("Documents","entradas","preci")

rasterPath = basepath.joinpath("raster")

vectorPath = basepath.joinpath("cadiz.shp")
*/

Creamos una variable denominada “dir_recorte” para crear el directorio en el que se van a almacenar los ráster recortados y establecemos una condición para que, si existe la carpeta “recortes” en ese directorio no haga nada, pero que si no existe una carpeta con ese nombre nos la genere automáticamente con el método .mkdir().

/**
dir_recorte = rasterPath.joinpath('recortes')

if not os.path.exists(dir_recorte):
    os.mkdir(dir_recorte)
*/

A continuación, creamos un bucle “for” para iterar la carpeta que contiene los ráster. Utilizamos el método “glob” para buscar aquellos ficheros con el formato .tif”. Dentro del bucle definimos la ruta de salida en cada momento de la iteración, dependiendo del ráster que se esté utilizando como entrada. El nombre que tomará cada ráster de salida será del tipo “Cadiz_[lyr-name]”, siendo “lyr.name” el nombre del ráster de entrada, incluyendo su extensión (.tif). Utilizamos el método .format(), el cual devuelve el argumento convertido en tipo texto.

Dentro del bucle también definimos los parámetros del algoritmo de “Clip raster by mask layer”. Hay que destacar que como ‘INPUT’ se va a utilizar el fichero ráster en cada momento de la iteración (variable “lyr”), que como ‘MASK’ se va a utilizar la capa vectorial, y que el ‘OUTPUT’ va a ser la ruta definida en el “output_path”. Añadimos al final de estas variables el método .as_posix(), el cual representa la cadena utilizando barras invertidas (/).

Tras los parámetros creamos la variable que almacenará el resultado, donde llamaremos al algoritmo con su ID («gdal:cliprasterbymasklayer») y con los parámetros como segundo argumento. Por último, añadiremos las capas resultantes al proyecto con QgsProject.instance().addMapLayer(). El nombre de las capas en el proyecto también estarán condicionadas por el bucle, por lo que el nombre también cambiará en función del ráster de entrada. Esto se puede generar utilizando .stem, el cual devuelve el nombre del fichero (en este caso el ráster de entrada) sin la extensión.

/**
for lyr in rasterPath.glob("*.tif"):
    output_path = dir_recorte.joinpath("Cadiz_{}".format(lyr.name))
    parameters = {
        'ALPHA_BAND': False,
        'CROP_TO_CUTLINE': True,
        'DATA_TYPE': 0,
        'INPUT': lyr.as_posix(),
        'KEEP_RESOLUTION': False,
        'MASK': vectorPath.as_posix(),
        'MULTITHREADING': False,
        'OPTIONS': '',
        'OUTPUT': output_path.as_posix(),
        'SET_RESOLUTION': False,
        'SOURCE_CRS': None,
        'TARGET_CRS': None,
        'X_RESOLUTION': None,
        'Y_RESOLUTION': None
    }

    resultado = processing.run("gdal:cliprasterbymasklayer", parameters)
    QgsProject.instance().addMapLayer(QgsRasterLayer(resultado['OUTPUT'], output_path.stem, 'gdal'))
*/

Como resultado del ejercicio tendremos en nuestra carpeta “recortes” los 10 ráster de precipitación (2010 a 2019) de la provincia de Cádiz.

Os comparto el código completo.

/**
from osgeo import gdal
import os
from pathlib import Path

basepath = Path().home().joinpath("Documents","entradas","preci")

rasterPath = basepath.joinpath("raster")

vectorPath = basepath.joinpath("cadiz.shp")

dir_recorte = rasterPath.joinpath('recortes')

if not os.path.exists(dir_recorte):
os.mkdir(dir_recorte)

for lyr in rasterPath.glob("*.tif"):
output_path = dir_recorte.joinpath("Cadiz_{}".format(lyr.name))
parameters = {
'ALPHA_BAND': False,
'CROP_TO_CUTLINE': True,
'DATA_TYPE': 0,
'INPUT': lyr.as_posix(),
'KEEP_RESOLUTION': False,
'MASK': vectorPath.as_posix(),
'MULTITHREADING': False,
'OPTIONS': '',
'OUTPUT': output_path.as_posix(),
'SET_RESOLUTION': False,
'SOURCE_CRS': None,
'TARGET_CRS': None,
'X_RESOLUTION': None,
'Y_RESOLUTION': None
}

resultado = processing.run("gdal:cliprasterbymasklayer", parameters)
QgsProject.instance().addMapLayer(QgsRasterLayer(resultado['OUTPUT'], output_path.stem, 'gdal'))
*/

(5 votos, promedio: 5,00 de 5)

Cargando...

Formación de calidad impartida por profesionales