The new version of QGIS is QGIS3 and it runs with Python 3 which introduces some change on the interaction with webservers with package “requests”. For those that are new to the IMERG images, those are some kind of the new TRMM images with precipitation estimation from multiple passive microwave (PMW) sensors on various precipitation-relevant satellites starting in March 2014. The IMERG images have a pixel resolution of 0.1 degrees and a temporal scale of 30 minutes; on the current panorama of precipitation estimates based on satellite-gauge, the IMERG data product with the highest spatial and temporal resolution available over the last 4 years.

This tutorial use Python to direct download IMERG precipitation images on a QGIS3 session. The tutorial is divided on three parts:

Download images from Nasa GesDisc
Conversion from HDF5 formato to Tiff
Add Tiff to QGIS canvas

Total amount of code used in this tutorial is no more than 60 lines on the three parts; Python code was done entirely with the standard QGIS3 and Python3 libraries. No external package was installed in order to download and translate the precipitation images.

Data manipulation

Nasa delivers IMERG precipitation images on HDF5 format that is not compatible with QGIS3, however these files could be opened with the GDAL library on the QGIS3 Python console. From the HDF5 layers, the precipitationCal dataset was the selected since it is recommended for general use by the IMERG developers. For the spatial analysis of large datasets, precipitation values would preferable to be manipulated as Numpy arrays rather than spatial raster, since the matrix algebra on Numpy allow us to make faster and easier additions, subtractions and filters on a large set of images as arrays.

See more technical details of IMERG precipitation estimates here:

https://pmm.nasa.gov/sites/default/files/document_files/IMERG_doc.pdf

Tutorial

First part: Data download

This part obtain the download links on the NASA Earthdata Search server (https://search.earthdata.nasa.gov/search) and store it on a local folder on the desktop. You need a Nasa Account that you can get for free on this link:

https://urs.earthdata.nasa.gov//users/new

A script in Python is run to download the links. Like other data products, the Imerg HDF name can have information about the dataproduct, algorithm version, date and provider, but this information could be useless for the scope of our analysis; this code also registers IMERG images with a particular file name pattern.

Since the urls are imported from a text file there are some breakline characters (‘\n’) that block the HDF request to the server. In order to get rid of these characters, the regular expression module (regex module or re) has been used in the tutorial to check each link for the presence of the ‘\n’ character and replace it with an empty space.

Video

Python code

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)

Second part: Raster format conversion

The second part is more intensive in Python commands and raster specification since it extracts one layer from each HDF file and store it as a separate Tiff file.

In order to make the matrix of precipitation values spatially adjusted we need to transpose the precipitation arrays, invert row order and set a geotransformation.

Note that there is a reassignment of the ‘raster’ value to None, this is relevant to close the file on the system. If we don't do this reassignment, the last raster on the loop will remain open and won't be available to be added to the canvas unless you exit QGIS.

Geotransformation resolution and extend depend on the raster metadata, for the tutorial the raster extension is the entire globe and the correspondent coordinates and resolution were previously setup on the geotransformation. For your analysis, you have to extract the raster resolution and specified coordinates of the raster dataset.

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

Third part: Add rasters to QGIS canvas

Video

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