La nueva versión de QGIS es QGIS 3 y se ejecuta con Python 3 que introduce algunos cambios en la interacción con los servidores web con paquetes de "requests". Para aquellos que son nuevos en imágenes IMERG, estas son algo así como las nuevas imágenes TRMM con estimación de precipitación de múltiples sensores de microondas pasivas (PMW) en varios satélites de precipitación relevantes cuya data comienza en Marzo del 2014. La imágenes IMERG tienen una resolución de pixel de 0.1 grados y una escala temporal de 30 minutos. En el actual panorama de las estimaciones de precipitación basadas en satélites, el producto de datos IMERG son las imágenes de más alta resolucion temporal y espacial disponibles de los últimos 4 años.
En este tutorial se usa Python para descargar directamente las imágenes de precipitación IMERG en una sesión de QGIS 3. El tutorial está dividido en tres partes:
- Descargar imágenes de la Nasa GesDisc.
- Conversión del formato HDF5 A TIFF.
- Añadir TIFF a la interfaz QGIS.
La cantidad total de código usado en este tutorial no es más de 60 líneas en las tres partes; el código Python se realizó con las librerías estándar de QGIS 3 y Python 3. No se instaló ningún paquete externo para descargar y convertir las imágenes de precipitación.
Manejo de datos
La Nasa ofrece imágenes de precipitacion IMERG en el formato HDF5 que no es compatible con QGIS 3, sin embargo, estos archivos se pueden abrir con las librerias GDAL en la consola Python QGIS 3. A partir de las capas HDF5, se seleccionó el dataset PrecipitacionCal ya que los desarrolladores de IMERG lo recomiendan para uso general. Para el análisis espacial de un gran conjunto de datos, sería preferible manejar los valores de precipitación como matrices Numpy en lugar de rasters espaciales, ya que la álgebra matricial en Numpy nos permite hacer sumas, restas y filtros en un gran conjunto de imágenes como matrices.
Ver más detalles técnicos de las estimaciones IMERG en este link:
https://pmm.nasa.gov/sites/default/files/document_files/IMERG_doc.pdf
Tutorial
Primera parte: Descarga de datos
En este parte se obtiene los enlaces de descarga en el servidor NASA Earthdata Search (https://search.earthdata.nasa.gov/search) y los guarda en una carpeta del escritorio. El usuario necesita tener una cuenta en Nasa, que puede obtener de forma gratuita en este enlace:
https://urs.earthdata.nasa.gov//users/new
Se ejecuta una secuencia de comandos en Python para descargar los enlaces. Al igual que otros productos, el nombre Imerg HDF puede tener información acerca de los datos del producto, versión del algoritmo, la fecha y el proveedor, pero esta información podría ser inútil para el alcance de nuestro análisis; este código también registra imágenes IMERG con un patrón particular de nombre de archivo.
Dado que las URL se importan desde una archivo de texto, hay algunos caracteres de quiebre de linea ('\ n') bloquean la solicitud de descarga al servidor. Para deshacerse de estos caracteres, se ha utilizado el módulo de expresión regular (módulo regex o re) en el tutorial para verificar la presencia del carácter '\ n' en cada enlace y reemplazarlo por un espacio vacío.
Video
Código de Python
import requests, os,re
os.chdir('C:\\Users\\USER\\Documents\\Infohataris\\DirectDownloadQGIS3\\Txt')
hdflinks = open('download.txt').readlines()
#clean links for breakline characters
for link in hdflinks:
index = hdflinks.index(link)
link = re.sub('\n','',link)
hdflinks[index] = link
#print(hdflinks)
for link in hdflinks:
print(link)
with requests.Session() as session:
req = session.request('get', link)
r = session.get(req.url, auth=('user', 'password'))
imagename = link.split('.')
with open('../HDF5/'+imagename[9][:14]+'.HDF5', 'wb') as f:
f.write(r.content)
Segunda parte: conversión de formato ráster
La segunda parte es más intensiva en los comandos de Python y la especificación de ráster, ya que extrae una capa de cada archivo HDF y la almacena como un archivo Tiff independiente.
Para hacer que la matriz de los valores de precipitación se ajuste espacialmente, tenemos que transponer las matrices de precipitación, invertir el orden de las filas y establecer una geotransformación.
Tenga en cuenta que hay una reasignación del valor 'ráster' a 'None', esto es relevante para cerrar el archivo en el sistema. Si no hacemos esta reasignación, el último ráster en el bucle permanecerá abierto y no estará disponible para ser agregado al lienzo a menos que salga de QGIS.
La resolución y extensión de geotransformación dependen de los metadatos del ráster, para el tutorial, la extensión de ráster es todo el globo y las coordenadas / resolución correspondientes se configuraron previamente en la geotransformación. Para su análisis, se debe extraer la resolución del ráster y las coordenadas especificadas de los metadatos del ráster.
Video
Code
import gdal, os, osr
from qgis.core import *
import qgis.utils
os.chdir('C:\\Users\\USER\\Documents\\Infohataris\\DirectDownloadQGIS3\\HDF5')
print(os.listdir(os.getcwd()))
totallinks = os.listdir(os.getcwd())
hdflinks = []
for link in totallinks:
if link[-4:] == 'HDF5':
hdflinks.append(link)
for link in hdflinks:
if link[-4:] == 'HDF5':
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
geotransform = ([-180,0.1,0,90,0,-0.1])
raster = gdal.GetDriverByName('GTiff').Create("../Tiff/"+link[:-5]+".tif",3600,1800,1,gdal.GDT_Float32)
raster.SetGeoTransform(geotransform)
srs=osr.SpatialReference()
srs.ImportFromEPSG(4326)
raster.SetProjection(srs.ExportToWkt())
raster.GetRasterBand(1).WriteArray(band_array.T[::-1])
raster = None
Tercera parte: Agregar rásters al lienzo QGIS
Vídeo
Code
import gdal, os, osr
from qgis.core import *
import qgis.utils
os.chdir('C:\\Users\\USER\\Documents\\Infohataris\\DirectDownloadQGIS3\\Tiff')
tiflinks = os.listdir(os.getcwd())
for link in tiflinks:
if link[-3:] == 'tif':
iface.addRasterLayer(link,link[:-4],"gdal")
