⚠️ This is the first operational version of the handbook, but it is still a work in progress and will be heavily updated during 2024! ⚠️

Hazard assessment for population exposed to drought#

Binder Droughts

Hazard assessment methodology#

In this workflow, we will assess the drought hazard data, expressed as Combined Drought Indicator (CDI).

Hazard data

The Combined Drought Indicator (CDI) that is implemented in the European Drought Observatory (EDO) is used to identify areas affected by agricultural drought, and areas with the potential to be affected. CDI can be downloaded from the Copernicus data server and is derived by combining three drought indicators produced operationally in the EDO framework - namely the

  • Standardized Precipitation Index (SPI),

  • the Soil Moisture Anomaly (SMA), and

  • the FAPAR Anomaly - in such a way that areas are classified according to three primary drought classes:

    1. “Watch” - indicating that precipitation is less than normal;

    2. “Warning” - indicating that soil moisture is in deficit; and

    3. “Alert” - indicating that vegetation shows signs of stress.

    Two additional classes - namely
    1. “Partial recovery” and

    2. “Recovery” - identify the stages of the vegetation recovery process.

Table 8 Classification scheme#

LEVEL

COLOUR

CLASIFICATION CONDITION

Watch

yellow

SPI-3 <-1 or SPI-1 < -2

Warning

orange

SMA <-1 and (SPI-3 < -1 or SPI-1 < -2)

Alert

red

\(\Delta\)FAPAR < -1 and (SPI-3 < -1 or SPI-1 < -2)

Partial recovery

brown

(\(\Delta\)FAPAR < -1 and (SPI-3m-1 < -1 or SPI-3 > -1) or
(\(\Delta\)FAPAR < -1 and (SPI-1m-1 < -2 or SPI-1 > -2)

Full recovery

grey

(SPI-3m-1 < -1 and SPI-3 > -1) or
(SPI-1m-1 < -2 or SPI-1 > -2)

Preparation#

Load libraries#

import os
import pooch

import numpy as np
import rasterio
from pathlib import Path
import rioxarray as rxr
import xarray as xr

import pyproj

import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

Create the directory structure#

In order for this workflow to work even if you download and use just this notebook, we need to set up the directory structure for the in- and output data.

Create a directory called drought_workflow in the same directory where this notebook is saved:

workflow_folder = 'simple_drought_workflow'
if not os.path.exists(workflow_folder):
    os.makedirs(workflow_folder)

Create a data directory inside the drought_workflow folder:

data_dir = os.path.join(workflow_folder,'data')
if not os.path.exists(data_dir):
    os.makedirs(data_dir)

Download the CDI#

In this workflow we will use data from the Copernicus Drought Observatory. The available datasets can be explored interactively, e.g., with the European Drought Observatory map viewer. Data can be obtained manually via the map viewer website (“Downloads” tab in the top right) or directly from their data directories. Here, we choose the latter and, as an example, download the CDI for the year 2022 as a NetCDF file into our data_dir folder:

filename = 'cdinx_m_euu_20220101_20221221_t.nc'
url = f'https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/DROUGHTOBS/Drought_Observatories_datasets/EDO_Combined_Drought_Indicator/ver3-0-1/{filename}'
known_hash = '7a40b16e5a0cfd53a22ef092f6a8bb511925972c9aaa9913a563c7e8585492e7'

pooch.retrieve(url=url, known_hash=known_hash, path=data_dir, fname=filename)

We can list all files in the data_dir directory using the os library.

with os.scandir(data_dir) as entries:
    for entry in entries:
        print(entry.name)
GHS_POP_E2020_GLOBE_R2023A_4326_30ss_V1_0.tif
9d9b95a2ca01c5c5a39324e22af110ce-GHS_POP_E2020_GLOBE_R2023A_4326_30ss_V1_0.zip
GHS_POP_E2020_GLOBE_R2023A_4326_30ss_V1_0.tif.ovr
GHS_POP_GLOBE_R2023A_input_metadata.xlsx
cdinx_m_euu_20220101_20221221_t.nc
GHSL_Data_Package_2023.pdf

Explore the CDI#

All downloaded netCDF files are stored in our data_dir, with file names starting with cdinx_m_euu_, followed by the start and end date period, and ending with _t.nc.

We take a look into our 2022 file downloaded earlier:

cdi_nc = xr.open_dataset(f"{data_dir}/cdinx_m_euu_20220101_20221221_t.nc")
cdi_nc
<xarray.Dataset> Size: 137MB
Dimensions:  (time: 36, lat: 950, lon: 1000)
Coordinates:
  * lat      (lat) int64 8kB 5497500 5492500 5487500 ... 762500 757500 752500
  * lon      (lon) int64 8kB 2502500 2507500 2512500 ... 7487500 7492500 7497500
  * time     (time) datetime64[ns] 288B 2022-01-01 2022-01-11 ... 2022-12-21
Data variables:
    cdinx    (time, lat, lon) float32 137MB ...
    3035     float32 4B ...
Attributes: (12/27)
    Source_Software:            dbinterface.py, dbexport.py, netcdf_handling.py
    creator_name:               Carolina Arias Munoz
    Conventions:                CF-1.6
    _CoordSysBuilder:           ucar.nc2.dataset.conv.CF1Convention
    date_created:               2023-06-23
    01.title:                   Combined Drought Indicator (CDI), v.3.0.1
    ...                         ...
    17.factsheet_url:           https://edo.jrc.ec.europa.eu/documents/factsh...
    18.jrc_data_catalogue_url:  https://data.jrc.ec.europa.eu/dataset/afa8a5e...
    19.sample_url:              /images/map_examples/dbio_data_previews/cdinx...
    20.metadata_last_updated:   2023-04-05
    21.values_legend:           0: No drought, 1: Watch. Precipitation defici...
    22.version_notes:           Current version: Version 3.0.1 covers data fr...

There are 36 time stamps for 12 months in the year 2022. This is because this index is calculated from satellite data that take 10 days to collect the information for the whole globe. The index of each time step corresponds to the start date of the associated 10-day time interval.

Tip

Explore the file content

Feel free to explore the content and structure of the datasets. Note the coordinates, dimensions and attributes!

Process the data#

Select the time of evaluation#

To check which times are available, you can inspect the data cells above, or print the time series like in the next cell:

cdi_nc.time
<xarray.DataArray 'time' (time: 36)> Size: 288B
array(['2022-01-01T00:00:00.000000000', '2022-01-11T00:00:00.000000000',
       '2022-01-21T00:00:00.000000000', '2022-02-01T00:00:00.000000000',
       '2022-02-11T00:00:00.000000000', '2022-02-21T00:00:00.000000000',
       '2022-03-01T00:00:00.000000000', '2022-03-11T00:00:00.000000000',
       '2022-03-21T00:00:00.000000000', '2022-04-01T00:00:00.000000000',
       '2022-04-11T00:00:00.000000000', '2022-04-21T00:00:00.000000000',
       '2022-05-01T00:00:00.000000000', '2022-05-11T00:00:00.000000000',
       '2022-05-21T00:00:00.000000000', '2022-06-01T00:00:00.000000000',
       '2022-06-11T00:00:00.000000000', '2022-06-21T00:00:00.000000000',
       '2022-07-01T00:00:00.000000000', '2022-07-11T00:00:00.000000000',
       '2022-07-21T00:00:00.000000000', '2022-08-01T00:00:00.000000000',
       '2022-08-11T00:00:00.000000000', '2022-08-21T00:00:00.000000000',
       '2022-09-01T00:00:00.000000000', '2022-09-11T00:00:00.000000000',
       '2022-09-21T00:00:00.000000000', '2022-10-01T00:00:00.000000000',
       '2022-10-11T00:00:00.000000000', '2022-10-21T00:00:00.000000000',
       '2022-11-01T00:00:00.000000000', '2022-11-11T00:00:00.000000000',
       '2022-11-21T00:00:00.000000000', '2022-12-01T00:00:00.000000000',
       '2022-12-11T00:00:00.000000000', '2022-12-21T00:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 288B 2022-01-01 2022-01-11 ... 2022-12-21
Attributes:
    standard_name:  time

For now we select just one time stamp for the CDI data. We can use the .sel() method of xarray to select one time step:

cdi = cdi_nc.sel(time='2022-08-01T00:00:00.000000000')

Select the area of interest#

We define the coordinates of the area of interest. Based on these coordinates we will be able to clip the dataset for further processing, and eventually display hazard and damage maps for the selected area.

To easily define an area in terms of geographical coordinates, you can use the Bounding Box Tool to select a region and get the corresponding coordinates. Make sure to select ‘CSV’ in the lower left corner and copy the values in the brackets into the cell below.

In this workflow example, we will concentrate on the area of Catalonia. It is found roughly between 0°E and 3.4°E longitude and 40.5°N and 42.9°N latitude:

xmin=0
ymin=40.5
xmax=3.4
ymax=42.9

However, the CDI dataset has its spatial coordinates in meters, so we first need to transform our coordinates to know where to clip to our selected area. For this we use the pyproj library. Here we use the proj-string of the Lambert Azimuthal Equal Area (LAEA) projection. This information can be found in the attributes of the CDI datatset for the 3035 variable:

cdi_nc["3035"].proj4_params
'+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000+y_0=3210000 +ellps=GRS80 +units=m +no_defs'
xyproj = pyproj.Proj('+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs')
latlonproj = pyproj.Proj('epsg:4326')
transformer = pyproj.Transformer.from_proj(latlonproj, xyproj)
(xmin_xyproj,ymin_xyproj) = transformer.transform(ymin, xmin)
(xmax_xyproj,ymax_xyproj) = transformer.transform(ymax, xmax)

Here are our transformed coordinates of the boundary box:

xmin_xyproj,ymin_xyproj
(3472008.40051275, 1989781.4970224823)
xmax_xyproj,ymax_xyproj
(3781134.3733767485, 2222944.719518429)

Caution

One very sneaky thing to note here is that latitudes (y dimension) in the dataset are decreasing. This is very important as this means we need to have the minimum and maximum y value in this order when selecting the area.

We can select our area in two ways: with the .sel() method of xarray or .rio.clip_box() provided by rasterio:

catalonia_cdi = cdi.sel(lat=slice(ymax_xyproj,ymin_xyproj), 
                        lon=slice(xmin_xyproj,xmax_xyproj))
catalonia_cdi = cdi.rio.clip_box(
    minx=xmin,
    miny=ymin,
    maxx=xmax,
    maxy=ymax,
   crs="EPSG:3050",
)

Let’s use xarray:

catalonia_cdi = cdi.sel(lat=slice(ymax_xyproj,ymin_xyproj), 
                        lon=slice(xmin_xyproj,xmax_xyproj))
catalonia_cdi
<xarray.Dataset> Size: 13kB
Dimensions:  (lat: 47, lon: 62)
Coordinates:
  * lat      (lat) int64 376B 2222500 2217500 2212500 ... 1997500 1992500
  * lon      (lon) int64 496B 3472500 3477500 3482500 ... 3772500 3777500
    time     datetime64[ns] 8B 2022-08-01
Data variables:
    cdinx    (lat, lon) float32 12kB ...
    3035     float32 4B ...
Attributes: (12/27)
    Source_Software:            dbinterface.py, dbexport.py, netcdf_handling.py
    creator_name:               Carolina Arias Munoz
    Conventions:                CF-1.6
    _CoordSysBuilder:           ucar.nc2.dataset.conv.CF1Convention
    date_created:               2023-06-23
    01.title:                   Combined Drought Indicator (CDI), v.3.0.1
    ...                         ...
    17.factsheet_url:           https://edo.jrc.ec.europa.eu/documents/factsh...
    18.jrc_data_catalogue_url:  https://data.jrc.ec.europa.eu/dataset/afa8a5e...
    19.sample_url:              /images/map_examples/dbio_data_previews/cdinx...
    20.metadata_last_updated:   2023-04-05
    21.values_legend:           0: No drought, 1: Watch. Precipitation defici...
    22.version_notes:           Current version: Version 3.0.1 covers data fr...

Reproject#

The only thing left to do before plotting is to reproject the CDI data. In order for rioxarray to be able to reproject the data, it needs to know the projection. NetCDF files store projection information in their variable or general attributes, so for rioxarray we need to write it using .rio.write_crs().

catalonia_cdi.rio.write_crs(3035, inplace=True)
catalonia_cdi_latlon = catalonia_cdi.rio.reproject("EPSG:4326")

Plot the CDI#

Basic plot#

With just one line code, we can create a simple plot of the dataset:

catalonia_cdi_latlon.cdinx.plot(cmap=plt.cm.PuRd, extend='max');
../../../../_images/f78021e62002b200d64accbaf14d3dc6e04b003442ac662607ddb5800db48915.png

For the plot to really help us to determine drought levels and in our climate risk assessement, we have to make a few changes to the map.

Custom plot#

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Set visible map section
ax.set_extent([xmin, xmax, ymin, ymax], crs=ccrs.PlateCarree())

# Add a land/ocean image in the background
ax.stock_img()

# Plot CDI data
cmap = mpl.colors.ListedColormap(["#ffffff", "#f0e442", "#f5aa00", "#dc050c","#0072b2","#cc79a7","#009e73", "#c8c8c8"])
c = catalonia_cdi_latlon.cdinx.plot(ax=ax, cmap=cmap, vmin=-0.5, vmax=7.5, add_colorbar=False, add_labels=False)

# Customize gridlines
gridlines = c.axes.gridlines(draw_labels=True, dms=False, xlocs=[0, 2, 4], ylocs=[40, 42])
# Add map features
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.RIVERS)

ax.set_title('Combined Drought Indicator', loc = "left")

# Create the categorical colorbar for the CDI data
cax_c = fig.add_axes([
    ax.get_position().x0,
    ax.get_position().y0,
    ax.get_position().width,
    0.08
])
cbar_c = fig.colorbar(c, orientation="horizontal", cax=cax_c)
cbar_c.ax.get_xaxis().set_ticks([])
labels = [
    'No drought',
    'Watch:\n rainfall\n deficit',
    'Warning:\n soil moisture\n deficit',
    'Alert','Recovery',
    'Temporary\n soil moisture\n recovery',
    'Temporary\n fAPAR\n recovery',
    'No data'
]
for j, label in enumerate(labels):
    cbar_c.ax.text(j, .5, label, ha='center', va='center', size=8.5)

# Add a marker and label at the location of Barcelona
barcelona_lon, barcelona_lat = 2.1686, 41.3874
c.axes.text(
    barcelona_lon,
    barcelona_lat,
    '+ Barcelona',
    horizontalalignment='left',
    fontsize=12,
    color='#000032',
    transform=ccrs.PlateCarree()
);
../../../../_images/a67409b1580d9175f9a01f433c5e5514305322c2ba9341253e23004984f9159a.png

Conclusions#

In this workflow, we have seen how to explore, process and visualise the drought data.

Contributors#

Milana Vuckovic, ECMWF
Maurizio Mazzoleni, Vrije Universiteit Amsterdam

Appendix I - Customised plot explained step by step#

We will:

  • Add coastlines, rivers and country borders using cartopy

  • Set distinct colour levels for CDI

  • Remove the default colorbar and labels

  • Add custom colorbars and title

  • Add custom gridlines

  • And finally add a label for Barcelona

Background#

First we make a basic plot and set geographical extent to our domain. This way we can zoom in without having to manipulate the data.
We can also add some (low resolution) background.

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Set visible map section
ax.set_extent([xmin, xmax, ymin, ymax], crs=ccrs.PlateCarree())

# Add a land/ocean image in the background
ax.stock_img();
../../../../_images/f33fb74c3e9d98435e6854888ba5e2d8febf5acff01595f172bd704bab3e87c8.png

Adding the data#

Next we add the data. We customize the levels on both plots. We will also remove the default colorbar and labels, so that we can add our own in the next step.

# Plot CDI data
cmap = mpl.colors.ListedColormap(["#ffffff", "#f0e442", "#f5aa00", "#dc050c","#0072b2","#cc79a7","#009e73", "#c8c8c8"])
c = catalonia_cdi_latlon.cdinx.plot(ax=ax, cmap=cmap, vmin=-0.5, vmax=7.5, add_colorbar=False, add_labels=False)

fig # show the updated figure in the notebook
../../../../_images/8c1918eb9d62493c28ca413359f8bf639ce1bd856ffd9e813dd873065ef618e6.png

Adding the geographical features#

Next we add the gridlines, coastlines, rivers and country borders.

# Customize gridlines
gridlines = c.axes.gridlines(draw_labels=True, dms=False, xlocs=[0, 2, 4], ylocs=[40, 42])
# Add map features
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.RIVERS)

fig # show the updated figure in the notebook
../../../../_images/3f96fc31052fe089ef11ca0161da142e903e6d716167390e98beb4cdf0ee9ffe.png

Adding the title#

Then we add the title.

ax.set_title('Combined Drought Indicator', loc = "left")

fig # show the updated figure in the notebook
../../../../_images/9ed353269e80e4ac27e62d470f2caff4be6b95895cabf60a709bfb6dd8d0e425.png

Adding the colorbar#

It is a little bit tricky to add colorbars to where we want them.
cax_p and cax_c variable define the positions and size of the colorbars. They get the default position and size of the colorbars and make them half as long and next to each other.
In the end we also make the tick labels smaller.

# Create the categorical colorbar for the CDI data
cax_c = fig.add_axes([
    ax.get_position().x0,
    ax.get_position().y0,
    ax.get_position().width,
    0.08
])
cbar_c = fig.colorbar(c, orientation="horizontal", cax=cax_c)
cbar_c.ax.get_xaxis().set_ticks([])
labels = [
    'No drought',
    'Watch:\n rainfall\n deficit',
    'Warning:\n soil moisture\n deficit',
    'Alert','Recovery',
    'Temporary\n soil moisture\n recovery',
    'Temporary\n fAPAR\n recovery',
    'No data'
]
for j, label in enumerate(labels):
    cbar_c.ax.text(j, .5, label, ha='center', va='center', size=8.5)

fig # show the updated figure in the notebook
../../../../_images/6d2ecffc134dc45ee6fa5800b85685e609ddc71c3c12e08039b453ad3e645fc0.png

Adding the custom text label#

Finally we can add the label for Barcelona.

# Add a marker and label at the location of Barcelona
barcelona_lon, barcelona_lat = 2.1686, 41.3874
c.axes.text(
    barcelona_lon,
    barcelona_lat,
    '+ Barcelona',
    horizontalalignment='left',
    fontsize=12,
    color='#000032',
    transform=ccrs.PlateCarree()
)

fig # show the updated figure in the notebook
../../../../_images/a67409b1580d9175f9a01f433c5e5514305322c2ba9341253e23004984f9159a.png