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

Extreme precipitation: Examples on how to carry out a climate risk assessment#

Click Binder to launch this notebook on MyBinder.
Click Heavy Rainfall to go to this workflow’s GitHub repository.

Examples of how to develop a risk assessment#

In this section, you can check real-world examples on how to carry out a risk assessment following the guidelines presented in this workflow. Feel free to explore the code and modify it according to your specific needs.

Prepare your workspace#

Load libraries#

# Libraries to download data and manage files
import os
import requests
import zipfile
import glob
import osmnx

# Libraries for numerical computations, array manipulation and statistics.
import xarray as xr
import pandas as pd
import numpy as np

# Libraries to handle geospatial data
import rioxarray as rio
import geopandas as gpd
from shapely import box
import pyproj
import rasterio

# Libraries to plot maps, charts and tables
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import from_levels_and_colors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import contextily as ctx

# Choosing the matplotlib backend
%matplotlib inline

Example 1: Are rainfall warning levels ready for climate change? A case study in Catalonia, Spain#

The region of Catalonia is located in the northest part of Spain. Due to enhanced convective rainfall events located over the mountainous coastal basins, Catalonia is highly vulnerable to extreme rainfall events that are concentrated between September and May. These events are becoming more frequent, with some resulting in severe urban floods with numerous casualties and significant socio-economic damages. For instance, the storm Gloria in 2020 caused widespread flooding across Catalonia, leaving extensive damage, casualties and more than 500 evacuations. Several areas recorded more than 180 mm of rainfall in 24 hours during the storm, beating historical records all over the region.

To increase the preparedness and response of the region, the Meteorological Service of Catalonia (SMC) issues warnings for Dangerous Meteorological Situations due to rainfall using two critical thresholds:

  1. Low risk: 100 mm/ 24 hours

  2. High risk: 200 mm/ 24 hours

These have been defined based on regional studies in the Catalonia region. However, considering the increased frequency of extreme precipitation events, the impacts associated with these and the number of triggered warnings per year, the regional authorities are interested in exploring how these current official thresholds could vary due to the influence of climate change.

The objective of this risk assessment is to evaluate how these critical rainfall thresholds defined by the SMC vary in the context of climate scenarios. The results can provide officials with valuable insights into the frequency shifts of these extreme events and the magnitude of changes that could be expected. This information is key for informed decision-making and enhancing preparedness measures to mitigate the impacts of climate change.

Step 1: Select your analysis area#

For this risk assessment, two analyses will be made:

  1. A specific location: We will focus on the city of Blanes, located in the province of Girona. Blanes has experienced extreme precipitation events in the past, including the recent storm Gloria mentioned earlier

  2. Regional assessment: We will zoom out our analysis to cover the entire Catalonia region. This extended scope will provide insights into the overall impact of extreme precipitation across Catalonia

bbox = [0, 40.2, 3.3, 43.5]; areaname = 'Catalonia'

Setting your directory structure#

The next cell will create the extreme_precipitation_risk folder in the directory where the notebook is saved. You can change the absolute path defined in the variable workflow_dir to create it elsewhere.

# Define the path for the extreme_precipitation_risk folder.
workflow_dir = 'extreme_precipitation_risk'
# workflow_dir = 'flashflood_workflow'
# Create the directory checking if it already exists.
os.makedirs(workflow_dir, exist_ok=True)

Now, we will create a subfolder for general data as well as specific subfolders to save data for our area of interest, plots generated by the code and new datasets computed within the workflow.

# Define directories for general data.
general_data_dir = os.path.join(workflow_dir, 'general_data')
# Define specific directories for the selected area
data_dir = os.path.join(workflow_dir, f'data_{areaname}')
results_dir = os.path.join(workflow_dir, f'results_{areaname}')
plots_dir = os.path.join(workflow_dir, f'plots_{areaname}')

# Create the directories checking if they already exist.
os.makedirs(general_data_dir, exist_ok = True)
os.makedirs(data_dir, exist_ok = True)
os.makedirs(results_dir, exist_ok = True)
os.makedirs(plots_dir, exist_ok = True)

We also need to define the path of the extreme_precipitation_hazard folder in order to access the files previously downloaded and created in the Hazard Assessment notebook.

# Define the hazard directory
hazard_dir = 'extreme_precipitation_hazard'
# Define the subfolder's paths
hazard_general_data_dir = os.path.join(hazard_dir, 'general_data')
hazard_data_dir = os.path.join(hazard_dir, f'data_{areaname}')
hazard_results_dir = os.path.join(hazard_dir, f'results_{areaname}')

Step 2: Calculate or identify your critical impact-based rainfall thresholds#

As mentioned above, the two critical rainfall thresholds identified are:

  1. Low risk: 100 mm/ 24 hours

  2. High risk: 200 mm/ 24 hours Where the frequency/return period of these thresholds varies across the Catalonia region. The local data used for this analysis is publicly provided by the SMC here.

For the city of Blanes, the critical impact-based rainfall threshold has been defined as 100 mm/24 hours for a 5-year return period.

Step 3: Load the Hazard data#

To carry out our risk assessment, we will use the hazard dataset generated by executing the code in the Hazard Assessment Notebook (Catalonia example). Thus, we will define the DATASET_HAZARD_ASSESSMENT as True. Our analysis will compare the baseline (observations) to the period of 2041-2070. For more details on how the hazard dataset used in this risk assessment was generated, please refer to the Hazard Assessment. Additionally, you can find information on the pre-calculated hazard datasets in the Hazard section of this workflow.

# Define global variable to define the dataset to use (False to use pre-computed datasets)
DATASET_HAZARD_ASSESSMENT = True

# Define global options for hazard dataset (matching the ones in Hazard Assessment or within the pre-computed dataset options)
GCM = 'ichec-ec-earth'
RCP = 'rcp85'

# Select which time-frames and durations within the available you want to study (Check for options above) or choose the ones computed in the Hazard Assessment.
TIME_FRAMES = {
    'historical' :'1976-2005',
    RCP : '2041-2070' 
}

# Define durations to read/download
DURATIONS = [3, 24]

# Dictionary of hazard dataset path's (filled automatically with the information of above)
HAZARD_PATHS = {}
# (if needed) Download pre-calculated European datasets
if not DATASET_HAZARD_ASSESSMENT:
    # Download the zip file
    print('Downloading data...')
    url_data = "http://gebrada.upc.es/climaax/extreme_precipitation_data.zip"
    tmp_file = "extreme_precipitation_data.zip"
    response = requests.get(url_data)
    open(tmp_file, "wb").write(response.content)

    # Unzip file to the data directory 
    print(f'...extracting data into {general_data_dir}...')
    with zipfile.ZipFile(tmp_file, 'r') as zip_ref:
        zip_ref.extractall(general_data_dir)

    # Delete zip file
    print('...removing .zip file.')
    os.remove(tmp_file)

# Fill the HAZARD_PATHS dictionary and (if needed) slice the european datasets to the area of interest 
for dur in DURATIONS:
    HAZARD_PATHS[dur] = {}
    for run, years in TIME_FRAMES.items():
        if DATASET_HAZARD_ASSESSMENT:
            HAZARD_PATHS[dur][run] = os.path.join(hazard_results_dir, f'idf_{dur}h_{GCM}_{run}_{years}.nc')

        else:
            path = glob.glob(os.path.join(general_data_dir, GCM, run, f'idf_{dur}h*_{years}.nc'))[0]
            path_out = os.path.join(data_dir, os.path.basename(path))
            if not os.path.exists(path_out):
                with xr.open_dataset(path, decode_cf = True) as ds:
                    # Cut bounding box area
                    # Reproject usual coordinates (in WGS84) to EURO-CORDEX's CRS.
                    CRS = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
                    transformer = pyproj.Transformer.from_crs('epsg:4326',CRS)
                    West, South, East, North = bbox
                    rlon_min, rlat_min, rlon_max, rlat_max = transformer.transform_bounds(South, West, North, East)
                    # slice European dataset to area of region using boundig box.
                    ds = ds.sel(rlat = slice(rlat_min, rlat_max), rlon = slice(rlon_min, rlon_max))
                    # Save in data_dir directory for further use
                    ds.to_netcdf(path_out)
            HAZARD_PATHS[dur][run] = path_out

Step 4: Specify your critical impact-based rainfall thresholds#

To carry our risk assessment, let’s specify the critical impact-based rainfall thresholds for both Blanes and Catalonia. We will define these in terms of Magnitude, Duration and Frequency.

Keep in mind: The code below helps us translate and link our real-world data to the climate dataset projections we’re using in this workflow (non-biased corrected EURO-CORDEX). This step is necessary for improving the robustness of our analysis with the climate projections.

def magnitude_frequency_shift(thresh_param, ds_baseline, ds_horizon):
    """ Returns the new magnitudes and frequencies of given threshold under climate change.

        thresh_param : list with the 3 elements of critical threshold [intensity, duration, frequency]
        ds_baseline : dataarray of expected precipitation of the baseline period. DataArray index is expected to be
                        just 'return_period'.
        ds_horizon : same specifications as ds_baseline for future horizon (e.g. 2070)
    """
    # Threshold parameters
    th_int, th_dur, th_freq = thresh_param

    # Access data (return_periods and values for baseline and future)
    ret_per_bs = ds_baseline.return_period.values
    ret_per_proj = ds_horizon.return_period.values
    prec_bs = ds_baseline.values
    prec_proj = ds_horizon.values

    # Translation of threshold into baseline eurocordex (interpolation of precipitation vs. log(return periods)).
    prec_base_aux = np.interp(np.log(th_freq), np.log(ret_per_bs), prec_bs)

    ### Change in frequency
    # Get new frequency by interpolation
    shift_ret_per = np.interp(prec_base_aux, prec_proj, np.log(ret_per_proj))
    shift_ret_per = np.ceil(np.exp(shift_ret_per))

    ### Change in magnitude using proxy of baseline.
    # Get new precipitation by interpolation
    prec_proj_aux = np.interp(np.log(th_freq), np.log(ret_per_proj), prec_proj)
    factor_prec = prec_proj_aux/prec_base_aux
    shift_prec = th_int * factor_prec

    # Rounding values to 2 decimals.
    shift_prec = round(shift_prec, 2)
    factor_prec = round(factor_prec, 2)

    return shift_ret_per, shift_prec, factor_prec

For a specific-site: Blanes city, Catalonia#

# Define the parameters of your threshold.
INTENSITY_TH = 100     # in mm
DURATION_TH = 24       # in hours
FREQUENCY_TH = 5       # in years

# Lat/lon coordinates for Blanes locality.
POINT = [41.674, 2.792]
# Transform POINT to Hazard dataset reference system
CRS = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
transformer = pyproj.Transformer.from_crs('epsg:4326', CRS)
rpoint = transformer.transform(POINT[0], POINT[1])

# Read historical IDF for point
ds_hist = xr.open_dataset(HAZARD_PATHS[DURATION_TH]['historical'], decode_cf=True)
ds_hist = ds_hist.sel(rlat = rpoint[1], rlon = rpoint[0], duration=DURATION_TH, method = 'nearest', drop=True).idf

# Read future IDF for point
ds_proj = xr.open_dataset(HAZARD_PATHS[DURATION_TH][RCP], decode_cf=True)
ds_proj = ds_proj.sel(rlat = rpoint[1], rlon = rpoint[0], duration=DURATION_TH, method='nearest', drop=True).idf

# Apply function to poing
shift_point = magnitude_frequency_shift([INTENSITY_TH, DURATION_TH, FREQUENCY_TH], ds_hist, ds_proj); del ds_hist, ds_proj

For the entire region of Catalonia#

print('Downloading Catalonia thresholds...')
path_rp_th = os.path.join(data_dir, 'RP_CAT_threshold_100mm.tiff')
url_data = "http://gebrada.upc.es/climaax/RP_CAT_threshold_100mm.tiff"
response = requests.get(url_data)
open(path_rp_th, "wb").write(response.content)

# Load and plot the return period for critical threshold.
print("Loading Catalonia thresholds...")
ds_rp_th = rio.open_rasterio(path_rp_th)

# Access projection and data values
proj_rp_th = ds_rp_th.rio.crs
transform_rp_rh = ds_rp_th.rio.transform()
nvar, ny, nx = ds_rp_th.shape
data = ds_rp_th.values

# Create new arrays for return period shift, precipitation shift and precipitation factor.
sh_rp = np.full((ny,nx), fill_value=ds_rp_th._FillValue, dtype='float32')
sh_prec = np.full((ny,nx), fill_value=ds_rp_th._FillValue, dtype='float32')
sh_prec_factor = np.full((ny, nx), fill_value=ds_rp_th._FillValue, dtype='float32')

# Transform lat/lon coordinates into Hazard Dataset CRS.
xy_th = np.asarray(np.where(data > 0))[[2,1],:]
trans_th_haz = pyproj.Transformer.from_crs(proj_rp_th, CRS)
lonlat_th = np.asarray(trans_th_haz.transform(ds_rp_th.x[xy_th[0]], ds_rp_th.y[xy_th[1]]))

# Read historical and future datasets
ds_hist = xr.open_dataset(HAZARD_PATHS[DURATION_TH]['historical'], decode_cf=True)
ds_future = xr.open_dataset(HAZARD_PATHS[DURATION_TH][RCP], decode_cf=True)

# Loop over each pixel to compute transformation
for i in range(len(xy_th[0])):
    th_rp_i = ds_rp_th.isel(x = xy_th[0, i], y = xy_th[1, i]).values[0]
    data_hist = ds_hist.idf.sel(rlon = lonlat_th[0,i], rlat = lonlat_th[1,i], method='nearest').sel(duration=DURATION_TH)
    data_future = ds_future.idf.sel(rlon = lonlat_th[0,i], rlat = lonlat_th[1,i], method='nearest').sel(duration=DURATION_TH)

    # Calling shift function
    shift = magnitude_frequency_shift([INTENSITY_TH, DURATION_TH, th_rp_i], data_hist, data_future)
    sh_rp[xy_th[1, i], xy_th[0, i]] = shift[0]
    sh_prec[xy_th[1, i], xy_th[0, i]] = shift[1]
    sh_prec_factor[xy_th[1, i], xy_th[0, i]] = shift[2]
def write_geotiff(filename, data):
    path_tiff = os.path.join(results_dir, filename)
    with rasterio.open(
        path_tiff,
        'w',
        driver='GTiff',
        height=data.shape[0],
        width=data.shape[1],
        count=1,
        dtype=data.dtype,
        crs=proj_rp_th.to_wkt(),
        nodata=ds_rp_th._FillValue,
        transform=transform_rp_rh,
        tiled=True,
        copy_src_overviews=True,
        compress='DEFLATE',
        predictor=3
    ) as dst:
        dst.write(data, 1)
        
# Save each variable as GeoTIFF in results directory.
filename_rp = f"RP_shift_threshold_100mm_{TIME_FRAMES['historical']}_{TIME_FRAMES[RCP]}.tiff"
write_geotiff(filename_rp, sh_rp)

filename_prec = f"PREC_shift_threshold_100mm_{TIME_FRAMES['historical']}_{TIME_FRAMES[RCP]}.tiff"
write_geotiff(filename_prec, sh_prec)

filename_factor = f"PREC_factor_threshold_100mm_{TIME_FRAMES['historical']}_{TIME_FRAMES[RCP]}.tiff"
write_geotiff(filename_factor, sh_prec_factor)

Step 5: Let’s explore our results!#

As a quick reminder, the objective of this risk assessment is to explore how the current critical impact-based rainfall thresholds will vary in terms of frequency and magnitude. In other words, we want to understand how rainfall patterns might change in the future while either maintaining the current frequency or magnitude of our critical thresholds.

For a specific-site: Blanes city#

In the case of Blanes (a specific location), we can print the results using the code below:

print(f"For the period of {TIME_FRAMES[RCP]}, the critical impact-based rainfall threshold of {INTENSITY_TH} mm/{DURATION_TH} hours, associated with a {FREQUENCY_TH}-year return period:")
print(f"   * If we want to maintain the same frequency (return period), the magnitude will vary by {(shift[2]-1)*100:.2f} % from the current threshold")
print(f"   * If we want to maintain the same magnitude, the frequency (return period) will change from {FREQUENCY_TH} to {shift_point[0]} years.")
For the period of 2041-2070, the critical impact-based rainfall threshold of 100 mm/24 hours, associated with a 5-year return period:
   * If we want to maintain the same frequency (return period), the magnitude will vary by 19.00 % from the current threshold
   * If we want to maintain the same magnitude, the frequency (return period) will change from 5 to 3.0 years.

For the entire region of Catalonia#

We can create some visualisations to better understand the changes in Magnitude and Frequency of our critical thresholds for the period of 2041-2070.

# Convert projection so matplotlib understands it.
proj_rp_map = ccrs.epsg(proj_rp_th.to_epsg())
# Path of rasters to plot.
r1 = path_rp_th
r2 = glob.glob(os.path.join(results_dir, 'RP_shift*.tiff'))[0]
# Define colorbar.
cmap_ret_per = mpl.cm.RdYlGn_r
# bounds_ret_per = [1, 2, 3, 4, 5, 10, 25, 50, 100, 200]
bounds_ret_per = [1, 5, 10, 20, 35, 50, 75, 100, 150, 200]
norm_ret_per = mpl.colors.BoundaryNorm(bounds_ret_per, cmap_ret_per.N)
fig_rp, ax_rp = plt.subplots(1,2, figsize = (18,11), subplot_kw = {'projection': proj_rp_map}, layout = 'constrained')
for i, file in enumerate([r1, r2]):
    ds_i = rio.open_rasterio(file, masked =True)
    im = ds_i.plot(ax = ax_rp[i], add_colorbar=False, norm = norm_ret_per, cmap=cmap_ret_per, alpha = 0.7)
    ctx.add_basemap(ax=ax_rp[i], crs=proj_rp_map, source=ctx.providers.OpenStreetMap.Mapnik,
               attribution = False)

ax_rp[0].set_title('Current return period')
ax_rp[1].set_title(f'Expected return period under {RCP}')
cb_rp = fig_rp.colorbar(im, ax = ax_rp[:], location='bottom', shrink = 0.7)
cb_rp.set_label('RP (years)', fontsize = 14)
fig_rp.suptitle(f"Shift of frequency for 100mm/24h threshold between in {TIME_FRAMES[RCP]}",fontsize = 16)
plt.show()

fig_rp.savefig(os.path.join(plots_dir, f'RP_shift_threshold_{TIME_FRAMES[RCP]}.png'))
../../../_images/5fee87a28f716acf2b4cfcd6d595c656e024332a6818bb67d568a3548fa9a707.png
# Path of rasters to plot.
r4 = glob.glob(os.path.join(results_dir, 'PREC_factor*.tiff'))[0]

# Define colormap for absolute values and relative change
cmap_rel = mpl.cm.RdYlGn
relative_bounds = np.arange(-50, 55, 5)
norm_rel = mpl.colors.BoundaryNorm(relative_bounds, cmap_rel.N, extend='both')


fig, ax = plt.subplots(figsize = (17,10), subplot_kw = {'projection': proj_rp_map}, layout = 'constrained')

ds_rel = rio.open_rasterio(r4, masked =True)
ds_rel = (ds_rel-1)*100
im_rel = ds_rel.plot(ax = ax, add_colorbar=False, norm = norm_rel, cmap=cmap_rel, alpha = 0.9)
cb_rel = fig.colorbar(im_rel, ax = ax, location ='bottom', shrink = 0.5 )
cb_rel.set_label('Precipitation shift (%)')
ctx.add_basemap(ax=ax, crs=proj_rp_map, source=ctx.providers.OpenStreetMap.Mapnik,
               attribution = False)

ax.set_title(f'Relative change under {RCP}')
fig.suptitle(f'Shift of magnitude for 100mm/24h threshold with fixed frequency in {TIME_FRAMES[RCP]}',fontsize = 16)
plt.show()
fig.savefig(os.path.join(plots_dir, f'PREC_shift_threshold_{TIME_FRAMES[RCP]}.png'))
../../../_images/114df7f8fac85831ab3b3c6bad70b29129967f5b2de81605dce73e4a4c8ebe10.png

Authors#

The Center of Applied Research in Hydrometrology from the Universitat Politècnica de Catalunya (CRAHI-UPC). The members involved in the development of this workflow were the following:

Erika R. Meléndez-Landaverde
Víctor González
Marta Gabarró Solanas
Daniel Sempere-Torres