Observing major flood events with GLOFAS

Observing major flood events with GLOFAS#

About#

This tutorial will demonstrate how to investigate a specific flood event using data from the GLObal Flood Awareness Service (GLOFAS). For this example we look at the major flood event in Bangladesh in June 2022 and assess the excess river discharge in comparison to the long-term (10-year) discharge statistics.

GloFAS is a product and service of the Copernicus Emergency Management Service. The Joint Research Centre of the European Commission is the entrusted entity responsible for CEMS GloFAS in terms of management, technical implementation and evolution. ECMWF is the designated contractor to implement the operational functionalities of the CEMS-COMP (Hydrological Forecast Computational Centre).

The tutorial is structured as follows:

(Install and import the required python packages)
Download the data of interest from the C3S Climate Data Store (CDS);
Explore the data with xarray, apply mask using the auxilary data, and produce demonstrative map plots;
Compute and plot time-series of the June 2022 river discharge against the 10-year period;
Compute and display a map of the number of days in June 2022 that the river discharge exceeded the 90th percentile of the 10-year period.

The steps and methods here serve as demonstrative examples and could be applied to other events and datasets, e.g. fire and heatwave events.

logo

Run the tutorial via free cloud platforms:

Setup your environment#

Before we begin we must prepare our environment. This includes installing the Application Programming Interface (API) of the CDS, and importing the various python libraries that we will need.

Install CDS API#

To install the CDS API, run the following command. We use an exclamation mark to pass the command to the shell (not to the Python interpreter).

!pip -q install cdsapi

# Other requirements:
!pip -q install cfgrib ecmwflibs
!pip -q install xarray netcdf4
!pip -q install matplotlib
!pip -q install cartopy

Import libraries#

We will be working with data in NetCDF format. To best handle this data we will use libraries for working with multidimensional arrays, in particular Xarray. We will also need libraries for plotting and viewing data, in this case we will use Matplotlib and Cartopy.

import os

# CDS API
import cdsapi

# Libraries for working with multidimensional arrays
import numpy as np
import xarray as xr

# Libraries for plotting and visualising data
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# Disable warnings for data download via API
import urllib3 
urllib3.disable_warnings()
# Disable xarray runtime warnings
import warnings
warnings.simplefilter("ignore", category=RuntimeWarning)

Search and download the data#

Enter your CDS API key#

We will request data from the Climate Data Store (CDS) programmatically with the help of the CDS API.

If you have not set up your CDS API credentials with a ~/.cdsapirc file, it is possible to provide the credentials when initialising the cdsapi.Client. To do this we must define the two variables needed to establish a connection: URL and KEY . The URL for the cds api is common and you do not need to modify that field. The KEY is string of characters made up of your your personal User ID and CDS API key. To obtain these, first register or login to the CDS (https://cds-beta.climate.copernicus.eu/), then visit https://cds-beta.climate.copernicus.eu/how-to-api and copy the string of characters listed after “key:”. Replace the ######### below with this string.

NOTE: If you have set up your cdsapi key using a ~/.cdsapirc you do not need to add your key below.

if os.path.isfile("~/.cdsapirc"):
    cdsapi_kwargs = {}
else:
    URL = 'https://cds-beta.climate.copernicus.eu/api'
    KEY = '##################################'
    cdsapi_kwargs = {
        'url': URL,
        'key': KEY,
    }

Here we specify a data directory in which we will download our data and all output files that we will generate:

DATADIR = './data_dir'
os.makedirs(DATADIR, exist_ok=True)

Find the river discharge data we want#

The historical GLOFAS data can be found in CDS catalogue entry.

We are going to look at the Bangladesh flood in June 2022. For our analysis we are going to compare June 2022 against all the June data from the 10-year reference period: 2012 to 2021. This period has been chosen for demonstrative purposes, a thorough analysis would choose a much longer reference period.

For our study we are going to make the following selection on the download page:

System version: version 4.0
Hydrological model: LISFLOOD
Product type: Consolidated
Variable: River discharge in the last 24 hours’
Year: 2012 to 2022
Month: June
Day: All days
Area: [30, 85, 20, 95,] = A bounding box for Bangladesh
Format: Grib

At the end of the download form, select “Show API request”. This will reveal a block of code, which you can simply copy and paste into a cell of your Jupyter Notebook (see cells below).

Download data#

… having copied the API request into the cell below, running this will retrieve and download the data you requested into your local directory. However, before you run the cell below, the terms and conditions of this particular dataset need to have been accepted in the CDS. The option to view and accept these conditions is given at the end of the download form, just above the “Show API request” option.

download_file = f"{DATADIR}/glofas-2012_2022.grib"
if not os.path.isfile(download_file):
    c = cdsapi.Client(
        **cdsapi_kwargs
    )
    c.retrieve(
        'cems-glofas-historical',
        {
            'system_version': 'version_4_0',
            'hydrological_model': 'lisflood',
            'product_type': 'consolidated',
            'variable': 'river_discharge_in_the_last_24_hours',
            'hyear': [f"{year}" for year in range(2012, 2023)],
            'hmonth': 'june',
            'hday': [f"{day:02d}" for day in range(1,31)],
            'format': 'grib',
            'area': [30, 85, 20, 95,],
        },
    ).download(download_file)

Compare the June 2022 river discharge to the previous 10 years#

We want to see how the June 2022 floods compare to the historical period accross the entire domain. We chose a histroical period of 10 years for demonstrative purposes, a thorough analysis would use a longer reference period.

The first step is to calculate some climatological statistics to represent our 10 year period (2012 to 2022).

Below we calculate the mean river discharge over the geospatial domain, then calculate 10 year mean, maximum, minimum and the 10th and 90th percentiles of the gespatial mean for each day in June.

We then calculate the geospatial mean for each day in June in 2022.

# select the 2012-2021 data
glofas_historical_period = glofas_data_masked.sel(time=slice("2012","2021"))

glofas_dis24_domain_mean = glofas_historical_period.dis24.mean(["longitude", "latitude"])
# calculate the mean discharge over the whole domain, for each day in June
glofas_xy_climatology = glofas_dis24_domain_mean.groupby('time.day').mean("time")

# calculate the minimum and maximum discharge over the whole domain, for each day in June
glofas_xy_clim_min = glofas_dis24_domain_mean.groupby('time.day').min("time")
glofas_xy_clim_max = glofas_dis24_domain_mean.groupby('time.day').max("time")
# Calculate the 10th and 90th percentiles of discharge rate of the whole domain
glofas_xy_clim_quantiles = glofas_dis24_domain_mean.groupby('time.day').quantile(dim="time", q=[0.1,0.9])

# calculate the mean discharge over the whole domain for each day in June 2022
glofas_xy_2022 = glofas_data_masked.dis24.sel(time="2022")
glofas_xy_2022_mean = glofas_xy_2022.mean(['longitude','latitude']).groupby('time.day').mean("time")

Plot these various statistics using matplotlib plotting methods.

# Create a figure and axis
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot()

# Plot the minimum to maximum range
ax.fill_between(
    np.arange(1, 31), glofas_xy_clim_min.values, glofas_xy_clim_max.values,
    color="k", alpha=0.2, label='min/max percentile range'
)

# Plot the 10th to 90th percentile range
ax.fill_between(
    np.arange(1, 31),
    glofas_xy_clim_quantiles.isel(quantile=0).values,
    glofas_xy_clim_quantiles.isel(quantile=1).values,
    color="k", alpha=0.5, label='10th/90th percentile range'
)

# Plot the mean for the 10 year period
(line_base,) = plt.plot(
    np.arange(1, 31), glofas_xy_climatology.values, color="k", label="Reference period mean"
)

# Plot the region where 2022 is greater than the mean in blue
ax.fill_between(
    range(1, int(glofas_xy_2022_mean.shape[0]) + 1),
    glofas_xy_climatology.values[: int(glofas_xy_2022_mean.shape[0])],
    glofas_xy_2022_mean.values,
    color="b",
    alpha=0.4,
    where=glofas_xy_2022_mean.values >= glofas_xy_climatology.values[: int(glofas_xy_2022_mean.shape[0])],
    interpolate=True, label="2022 excess discharge"
)

# Plot the region where 2022 is less than the mean in red
ax.fill_between(
    range(1, int(glofas_xy_2022_mean.shape[0]) + 1),
    glofas_xy_climatology.values[: int(glofas_xy_2022_mean.shape[0])],
    glofas_xy_2022_mean.values,
    color="r",
    alpha=0.4,
    where=glofas_xy_2022_mean.values < glofas_xy_climatology.values[: int(glofas_xy_2022_mean.shape[0])],
    interpolate=True,
)

# Plot the 2022 data points
dots = ax.scatter(
    range(1, int(glofas_xy_2022_mean.shape[0]) + 1),
    glofas_xy_2022_mean.values,
    color="k",
    # label="2022",
    zorder=2,
)

# Add a legend and axis labels
ax.legend(loc='upper left')
plt.xlabel("Day in June")
plt.ylabel("Daily river discharge ($m^{3}s^{-1}$)")

Text(0, 0.5, 'Daily river discharge ($m^{3}s^{-1}$)')

_images/bc7a7f44d5527d69bf94f705642153ebeb50f0afb69da435676865cc78a8391f.png

From this plot, we see that around the 13th of June the 2022 floods began to exceed the 10-year historical maximum, and this continued until the 22nd.

However, this plot doesn’t give us an idea of how the severity of the flooding varied spatially, so let’s have a look at that.

Map the number of days that the river discharge exceeded the 90th percentile of the 10-year period.#

We can also look at the spatial distribution of this data and xarray provides some useful tools for indexing calcualating various statisitics. Below we calculate the 10th and 90th percentiles accross the domain then compute where the river discharge exceeded the 90th percentile. This is a boolean data array where 1 indicates the 90th percentile was exceeded and 0 indicates it was not.

There are a number of commented out commands which can be used to calculate other statistics which you can investigate at your leisure.

# Calculate the climatology for the historical period: 2012-2021
glofas_historical_period = glofas_data_masked.sel(time=slice("2012","2021"))

# glofas_climatology = glofas_data_masked.groupby('time.dayofyear').mean("time")
glofas_climat_quantiles = glofas_historical_period.groupby('time.dayofyear').quantile(dim="time", q=[0.1,0.9])

# # Calculate the anomaly w.r.t the above climatology
# glofas_anomaly = glofas_data_masked.groupby('time.dayofyear') - glofas_climatology
# # Calculate this a relative anomaly
# glofas_relanom = glofas_anomaly.groupby('time.dayofyear')/glofas_climatology

glofas_high = glofas_data_masked.groupby('time.dayofyear') > glofas_climat_quantiles.sel(quantile=0.9)
glofas_high = glofas_high.where(glofas_data.uparea>250e6)

We can then plot the map using our predefined plotting function

# Plot the data:
year='2022'
number_of_high_days = glofas_high.sel(time=year).dis24.sum(dim='time')
number_of_high_days = number_of_high_days.assign_attrs(glofas_high.dis24.attrs)
plot_map(
    number_of_high_days, vmax=25, cmap='Purples', cbar_label='Number of days'
)

_images/4abed668d137956ffa55fe9d6e2d41e536bd1207bd01123a41ca9864f1ba7cb0.png

We can see that some areas to the northeast of Dhaka experienced discharge exceeding the 90th percentile for over 20 days in June.

Observing major flood events with GLOFAS

Contents

Observing major flood events with GLOFAS#

About#

Setup your environment#

Install CDS API#

Import libraries#

Search and download the data#

Enter your CDS API key#

Find the river discharge data we want#

Download data#

Explore and view the data#

Getting an overview of the study area#

Accessing the upstream area data#

Reduce the upstream area to the domain of the data#

Mask the river discharge data#

Compare the June 2022 river discharge to the previous 10 years#

Map the number of days that the river discharge exceeded the 90th percentile of the 10-year period.#