Tutorial de Operaciones Masivas en Rásters con QGIS 3 y Python

Existen nuevas herramientas y recursos disponibles para comprender el cambio climático, la dinámica del uso del suelo, el ciclo del agua y otras partes de nuestro entorno físico. Muchos datos espaciales vienen en formato raster y están disponibles en servidores web, que tienen imágenes registradas cada año, mes, día, hora, media hora o minuto. Si queremos evaluar un fenómeno físico debemos analizar un gran conjunto de datos.

El enfoque y las herramientas SIG tradicionales se centraron en realizar un algoritmo espacial en una sola capa o en un pequeño grupo de capas, sin embargo, necesitamos desarrollar y usar herramientas para aumentar nuestras capacidades de manipulación de datos a un número indefinido de capas; afortunadamente esto se puede hacer con QGIS y la consola Python.

Si manejamos datos ráster es importante recordar que todos los rásteres son matrices con una georreferenciación espacial. Cuando la cantidad de operaciones requerida para los rásters es enorme, es preferible convertir los rásters en matrices (o arreglos 2D Numpy), realizar operaciones con herramientas Numpy y luego geotransformarlas en matrices. Con respecto a la idoneidad espacial de nuestra álgebra matricial, todos los archivos ráster deben tener la misma resolución, proyección y dimensión.

Código de Python

Python es capaz de hacer múltiples operaciones en imágenes rásters y no se pretende que los "scripts" sean largos, estando en el orden de 40 a 60 líneas dependiendo del nivel del análisis requerido. La mayor eficiencia de Python para el análisis de rasters proviene de los bucles, condicionales, operaciones en Numpy y herramientas para abrir y escribir archivos. Hay dos versiones del código de Python para este tutorial, una para archivos HDF5 y para archivos TIFF; ambas que funcionan casi en la misma dirección. Aquí está la versión HDF5:

import gdal, os, osr
import numpy as np
from qgis.core import *
import qgis.utils

os.chdir('C:\\Users\\USER\\Documents\\Infohataris\\MassiveOperationwithRasterswithPyQGIS\\HDF5')
print(os.listdir(os.getcwd()))
totallinks = os.listdir(os.getcwd())

hdflinks = []
for link in totallinks:
    if link[-4:] == 'HDF5':
        hdflinks.append(link)

sum_array = np.zeros([1800,3600])

for link in hdflinks:
    hdf_ds = gdal.Open(link, gdal.GA_ReadOnly)
    print(link)
    band_ds = gdal.Open(hdf_ds.GetSubDatasets()[5][0], gdal.GA_ReadOnly) #choose the 5th dataset that corresponds to precipitationCal
    band_array = band_ds.ReadAsArray()
    band_array[band_array<0] = 0 #filter all NaN values that appear as negative values, specially for the tiff representation       
    sum_array = band_array.T[::-1]*0.5 + sum_array

geotransform = ([-180,0.1,0,90,0,-0.1])
raster = gdal.GetDriverByName('GTiff').Create("../Sum/SumImagesHdf.tif",3600,1800,1,gdal.GDT_Float32)
raster.SetGeoTransform(geotransform)
srs=osr.SpatialReference()
srs.ImportFromEPSG(4326)
raster.SetProjection(srs.ExportToWkt())
raster.GetRasterBand(1).WriteArray(sum_array)
raster = None

iface.addRasterLayer("../Sum/SumImagesHdf.tif","SumImagesHdf.tif","gdal")

Tutorial

Este tutorial desarrolla la suma de 48 archivos raster con algunas operaciones matriciales como transposición y orden inverso de filas. El conjunto de datos proviene de las imágenes de precipitación IMERG cada media hora durante un día (1 de enero de 2017).

Datos de entrada

Descarga los datos de entrada aquí.