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