Using PyQGIS to convert GRIB files to GeoTiff and add to QGIS project

Introduction

Using the Python API can extend flexibility through automation when working with data in GIS. It allows users access to the entire Python ecosystem of modules to for analysis, the ability to create standalone scripts using PyQGIS, run custom processing routines from within a QGIS project, and even create custom plugins.

Working with scientific data formats can present challenges when it comes to using those data within GIS. The files may have multiple dimensions, or their georeferencing information may not be recorded in the file in a way that QGIS understands.

Using Python scripts written to automate the conversion of these file formats to GeoTiff, can be written once and used over and over again, which can save time in the long run.

In this tutorial, we will provide an example script for converting a GRIB file to Geotiff, and adding this to a QGIS project.

Learning outcomes

• Familiarity with how to run Python scripts from within a QGIS project.

• An understanding of the potential for using Python and QGIS together.

• Awareness of how you can use the Python Core library, as well as any other Python modules, together with QGIS commands to add flexibility to your processing routines.

Tools required

QGIS Desktop v3.x. We recommend using the current Long Term Release version.

Dataset

The GRIB file we will be using is ERA5 reanalysis, hourly sea surface temperature and surface pressure for the 13-17th of May 2023.

Reference: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview

We will also be using a coastline dataset of the world, accessed from NaturalEarth.com

You can download these datasets from the link below:

https://files.bas.ac.uk/magic_products/v1/1037--geospatial-training/12--data-bundles/qgis-python-bundle.zip

Unzip the coastline data, ie. ne_110m_coastline.zip.

The Exercise

Open a new QGIS Project, and load the ne_110m_coastline layer.

Basic command line

We will start by demonstrating how to access the Python prompt and use this as a command line.

In your project, open the Python prompt. This is accessible either via the icon on the toolbar, or via the Main Menu > Plugins > Python Console.

QGIS User Interface with Python Console activated. Red circle highlights toolbar icon to activate Console.

Type the following command in to the Python prompt:

print('Hello, World!')

As you will see, the QGIS Python prompt acts exactly as an interactive Python command line does. From here, you can type single commands and interact with layers in your project, and access the Python Core library, and other modules as required.

Next, type the following in to the prompt:

iface.mapCanvas().layers()

This command creates a Python list containing references to all the layers which form part of your QGIS project.

The PyQGIS cookbook is a great resource for exploring what you can do with the API: https://docs.qgis.org/3.28/en/docs/pyqgis_developer_cookbook/index.html

Running a standalone script

Next, we'll load in a script (GRIB_to_GTIFF.py) to run within your project. This will convert the ERA5_hourly_SST_SP_May_2023.grib file to GeoTiff format and demonstrate how to load the converted files in to your project.

First, locate the tool in the Python prompt called Show Editor. Next, load the script GRIB_to_GTIFF.py and edit the following line to reflect the folder path where you have saved ERA5_hourly_SST_SP_May_2023.grib.

os.chdir(r"C:/training/")

QGIS interface showing Python prompt and script editor

Have a read through the script and the code comments to get an understanding of what the script will do.

The GRIB file contains many layers, each representing an hour of data. The script will convert these to GeoTiff files, and then calculate a daily average for both sea surface temperature and surface pressure.

The final two lines of the script will add the daily averages to your QGIS project, ie:

iface.addRasterLayer("SSTK_2023-05-13.tif", "SSTK 2023-05-13")
iface.addRasterLayer("SP_2023-05-13.tif", "SP 2023-05-13")

Go ahead and run the script, using the Run Script tool in the toolbar of the Python prompt.

Running the script should result in two layers being added to your project. You will also see that there are a number of GeoTiff files added to your working directory. These are the individual GRIB layers, which have been converted to GeoTiff.

What next?

Have a go at editing the script to experiment with what else you can do. Some ideas to get you started:

• Amend the script to calculate a daily average for a different day

• Following the step where the script adds the daily averages to the project, add some code to change the styling of those files