The return period of extreme precipitation

In this tutorial, we’ll explore what a return period is and how we can use it to define extreme precipitation. We will use historical climate model outputs to calculate the thresholds of daily precipitation associated with return periods of, e.g., 10-, 20-, and 50-years. To do so, we will use preprocessed climate model outputs for the Molise region in Italy (NUTS2 region ITF2).

Settings¶

Modules used¶

import pandas as pd
#import matplotlib.pyplot as plt
#import seaborn as sns
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import numpy as np
import os
from pathlib import Path
from scipy.stats import genextreme

# Configure plotting
# plt.rcParams['figure.figsize'] = (12, 6)
# plt.rcParams['font.size'] = 11

User settings¶

# Configuration
admin_id = 'ITF2'

# Reference period for baseline statistics
baseline_start = 1976
baseline_end = 2005

# Analysis settings
rolling_window = 30  # years for rolling mean

# Paths
workdir = Path("/home/nejk/code/extreme_precip_exposure")
os.chdir(workdir)
data_dir = workdir / 'data' / admin_id / 'extreme_precip_exposure'

# Load preprocessed data / Sample datasets
print(f"Region: {admin_id} (Molise, Italy)")
sample_file = data_dir / f'{admin_id}_precip_hist_sample_model_gridpoint.csv'
print(f"Sample data: {sample_file}")
sample_model_file = data_dir / f'{admin_id}_precip_hist_amax_sample_model.csv'
print(f"Sample model data: {sample_model_file}")
sample_point_file = data_dir / f'{admin_id}_precip_hist_amax_sample_gridpoint.csv'
print(f"Sample point data: {sample_point_file}")

# Load preprocessed data / Full datasets
hist_file = data_dir / f'{admin_id}_precip_hist_amax_all_points.csv'
print(f"Historical data: {hist_file}")
proj_file = data_dir / f'{admin_id}_precip_proj_amax_all_points.csv'
print(f"Projection data: {proj_file}")

Region: ITF2 (Molise, Italy)
Sample data: /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_sample_model_gridpoint.csv
Sample model data: /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_amax_sample_model.csv
Sample point data: /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_amax_sample_gridpoint.csv
Historical data: /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_amax_all_points.csv
Projection data: /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_proj_amax_all_points.csv

Data used¶

For this tutorial, we’ll use preprocessed data, which we downloaded and preprocessed for the Molise region (Italy) in the

how-to guide on processing precipitation from climate model outputs.

In this how-to-guide, we downloaded the CORDEX-EUR11 ensemble from the Gridded datasets underpinning the Copernicus Interactive Climate Atlas and more specifically, the variable monthly maximum 1-day precipitation.

The data that we load here represents either one (random) model but all points in the region and timesteps available, loaded as test_file, or all models but only the annual and regional maximumm (stored in hist_file and proj_file for the historical and future periods, respectively).

Return periods explained¶

Example timeseries of monthly maximum of daily precipitation¶

We’ll use the outputs from one (random) model here to illustrate a timeseries of monthly maximum of daily precipitation.

# Read test_file and plot timeseries of monthly maximum precipitation for one random point in the region
point_df = pd.read_csv(sample_file, parse_dates=['time'])
# remove duplicate rows
point_df = point_df.drop_duplicates()
point_df['year'] = point_df['time'].dt.year
point_df.head()
print(f"The sample data contains {len(point_df)} rows and covers the period from {point_df['time'].min().date()} to {point_df['time'].max().date()}. ")
print(f"\nThe underlying climate model is entitled {point_df['member_id'].unique()[0]}.")
print(f"The naming indicates the global and regional climate model (GCM and RCM, respectively) institutions, models, and variants:")
print(f" - GCM institution: {point_df['gcm_institution'].unique()[0]}")
print(f" - GCM model: {point_df['gcm_model'].unique()[0]}")
print(f" - GCM variant: {point_df['gcm_variant'].unique()[0]}")
print(f" - RCM institution: {point_df['rcm_institution'].unique()[0]}")
print(f" - RCM model: {point_df['rcm_model'].unique()[0]}")
print(f" - RCM variant: {point_df['rcm_variant'].unique()[0]}")

# Count the number of unique grid points in the sample data using the 'lon' and 'lat' column combinations
unique_points = point_df[['lon', 'lat']].drop_duplicates()
num_unique_points = len(unique_points)
print(f"\nThe sample data contains {num_unique_points} unique grid point(s).")

The sample data contains 432 rows and covers the period from 1970-01-01 to 2005-12-01. 

The underlying climate model is entitled CCCma_CanESM2_r1i1p1_CLMcom_CCLM4-8-17_v1.
The naming indicates the global and regional climate model (GCM and RCM, respectively) institutions, models, and variants:
 - GCM institution: CCCma
 - GCM model: CanESM2
 - GCM variant: r1i1p1
 - RCM institution: CLMcom
 - RCM model: CCLM4-8-17
 - RCM variant: v1

The sample data contains 1 unique grid point(s).

fig = go.Figure()

fig.add_trace(go.Scatter(
    x=point_df['time'],
    y=point_df['rx1day'],
    mode='lines',
    name='Monthly max daily precipitation',
    line=dict(color='steelblue', width=1.5),
))

annual_maxima = point_df.copy()
annual_maxima['annual_max'] = annual_maxima.groupby('year')['rx1day'].transform('max')
annual_maxima = annual_maxima[annual_maxima['rx1day'] == annual_maxima['annual_max']]

fig.add_trace(go.Scatter(
    x=annual_maxima['time'],
    y=annual_maxima['annual_max'],
    mode='markers',
    name='Annual Maximum',
    marker=dict(size=7, color='steelblue', symbol='circle'),
))

fig.update_layout(
    title=dict(text=f'Monthly Maximum of Daily Precipitation<br><sup>{point_df["member_id"].unique()[0]},lon={point_df["lon"].unique()[0]}, lat={point_df["lat"].unique()[0]}</sup>'),
    xaxis=dict(title='Time', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    yaxis=dict(title='Precipitation [mm]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=0.6, y=1, xanchor='left', yanchor='top'),
    width=800,
    height=500,
)

fig.show()

Through the selection of the monthly maximum of daily precipitation, we have already performed a first screening for extremes at the daily scale. However, we can also perform a second screening and look at annual maxima. These points are already highlighted in the plot above...

Using annual maxima to define extreme events¶

Using this data, we can, for example, understand during which month the most intense daily precipitation falls or identify months during which precipitation extremes are rare.

point_df = point_df.copy()
point_df.loc[:, 'month'] = point_df['time'].dt.month
point_df.loc[:, 'year'] = point_df['time'].dt.year

# Print statistics for this point, indicating months with highest and lowest precipitation, and the overall mean and standard deviation
mean_precip = point_df['rx1day'].mean()
std_precip = point_df['rx1day'].std()
max_precip = point_df['rx1day'].max()
min_precip = point_df['rx1day'].min()
max_month = point_df.loc[point_df['rx1day'].idxmax(), 'month']
min_month = point_df.loc[point_df['rx1day'].idxmin(), 'month']
print(f"Statistics for point ({point_df['lon'].iloc[0]}, {point_df['lat'].iloc[0]}):")
print(f" - Mean monthly maximum daily precipitation: {mean_precip:.2f} mm")
print(f" - Standard deviation: {std_precip:.2f} mm")
print(f" - Maximum monthly maximum daily precipitation: {max_precip:.2f} mm (in month {max_month})")
print(f" - Minimum monthly maximum daily precipitation: {min_precip:.2f} mm (in month {min_month})")

Statistics for point (14.5625, 41.6875):
 - Mean monthly maximum daily precipitation: 14.62 mm
 - Standard deviation: 12.86 mm
 - Maximum monthly maximum daily precipitation: 80.00 mm (in month 9)
 - Minimum monthly maximum daily precipitation: 0.00 mm (in month 7)

# Annual cycle plot, highlight the annual maximum per year using a colored dot in the same color as the year lines
month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
annual_max_df = point_df.groupby('year').apply(lambda x: x.loc[x['rx1day'].idxmax()], include_groups=False).reset_index()
mean_cycle = point_df.groupby('month')['rx1day'].mean().reset_index()

years = sorted(point_df['year'].unique())
colors = [f'hsl({int(240 - 240 * (y - years[0]) / max(1, years[-1] - years[0]))}, 70%, 50%)' for y in years]
year_color = dict(zip(years, colors))

fig = go.Figure()

for year, color in zip(years, colors):
    year_data = point_df[point_df['year'] == year].sort_values('month')
    fig.add_trace(go.Scatter(
        x=year_data['month'],
        y=year_data['rx1day'],
        mode='lines',
        line=dict(color=color, width=1.25),
        opacity=0.5,
        name=str(year),
        showlegend=False,
        hovertemplate=f'{year}<br>Month: %{{x}}<br>%{{y:.1f}} mm<extra></extra>',
    ))

# Annual maxima dots
fig.add_trace(go.Scatter(
    x=annual_max_df['month'],
    y=annual_max_df['rx1day'],
    mode='markers',
    marker=dict(
        color=annual_max_df['year'],
        colorscale='Viridis',
        size=8,
        showscale=True,
        colorbar=dict(title='Year', thickness=12),
    ),
    name='Annual Maximum',
    hovertemplate='Year: %{marker.color}<br>Month: %{x}<br>%{y:.1f} mm<extra></extra>',
))

# Mean annual cycle
fig.add_trace(go.Scatter(
    x=mean_cycle['month'],
    y=mean_cycle['rx1day'],
    mode='lines',
    line=dict(color='black', width=3),
    name='Mean Annual Cycle',
))

fig.update_layout(
    title=dict(text=f'Annual Cycle of Monthly Maximum Daily Precipitation<br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].iloc[0]}, lat={point_df["lat"].iloc[0]}</sup>'),
    xaxis=dict(
        title='Month',
        tickmode='array',
        tickvals=list(range(1, 13)),
        ticktext=month_names,
        showgrid=True, gridcolor='lightgrey', gridwidth=0.5,
    ),
    yaxis=dict(title='Precipitation [mm]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=0.025, y=0.95, xanchor='left', yanchor='top'),
    width=800,
    height=500,
)

fig.show()

Can you spot the annual maxima for each year in this plot?

No single year captures all deviations from the mean — that’s why we need a long climatology to define an extreme event. The World Meteorological Organisation suggests to use (at least) 30 years of data to do so!

# print months with most and least annual maximum precipitation
annual_max_idx = point_df.groupby('year')['rx1day'].idxmax()
annual_max_months = point_df.loc[annual_max_idx].set_index('year')['month']
most_common_max_month = annual_max_months.mode().iloc[0]
# account for months without any annual maximum precipitation (i.e., if a month has no annual maximum precipitation, it should be counted and listed as least common month!)
all_months = pd.Series(np.arange(1, 13))
least_common_max_month = all_months[~all_months.isin(annual_max_months)].tolist()
if not least_common_max_month:
    least_common_max_month = annual_max_months.value_counts().idxmin()

print(f"The most common month for the annual maximum precipitation is month {most_common_max_month} with {annual_max_months.value_counts().max()} occurrences.")  
print(f"The least common month for the annual maximum precipitation are months {least_common_max_month}.")

The most common month for the annual maximum precipitation is month 10 with 9 occurrences.
The least common month for the annual maximum precipitation are months [4, 7, 8].

A statistical model for the occurrence of extremes¶

The annual maxima help us to further refine what we call an extreme precipitation event — and they are the basis for the definition of return periods, which are expressed in units of years.

So, let’s look at these annual maxima a bit more closely, and let’s plot them on a histogram to understand their occurrence. To model the behavior of the annual maxima of daily precipitation at this grid point, we also fit a Generalised Extreme Value distribution to the histogram.
The GEV has three parameters that determine the tail behaviour of the distribution: location ( $\mu$ ), scale ( $\sigma$ ), and shape ( $\xi$ ).

# Fit a GEV distribution to the annual maxima and plot it on top of the above histogram
gev_params = genextreme.fit(annual_max_df['rx1day'])
x = np.linspace(annual_max_df['rx1day'].min(), annual_max_df['rx1day'].max(), 100)
pdf = genextreme.pdf(x, *gev_params)

# Build histogram counts manually so Plotly can normalise to density
counts, bin_edges = np.histogram(annual_max_df['rx1day'], bins=20, density=True)
bin_width = bin_edges[1] - bin_edges[0]

fig = go.Figure()

fig.add_trace(go.Bar(
    x=(bin_edges[:-1] + bin_edges[1:]) / 2,
    y=counts,
    width=bin_width * 0.95,
    marker=dict(color='steelblue', line=dict(width=0)),
    name='Annual max precipitation',
))

fig.add_trace(go.Scatter(
    x=x,
    y=pdf,
    mode='lines',
    line=dict(color='indigo', width=2),
    name='Fitted GEV PDF',
))

fig.update_layout(
    title=dict(text=f'Histogram of Annual Maximum of Daily Precipitation with Fitted GEV<br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].iloc[0]}, lat={point_df["lat"].iloc[0]}</sup>'),
    xaxis=dict(title='Precipitation [mm]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    yaxis=dict(title='Density', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=1, y=1, xanchor='right', yanchor='top'),
    bargap=0,
    width=800,
    height=500,
)

fig.show()

Definition of return periods¶

The above histogram and fitted GEV distribution already indicate probability of an event to occur, e.g., an event larger than 80 mm per day (intensity), occurring. The return period is defined using this exceedance probability.

Using the fitted GEV, we can use its probability density function (PDF) and cumulative density function (CDF) to estimate the probability of exceedance for specific thresholds:

The return period $T$ of an extreme precipitation event is defined as the average time between two events of the same or greater intensity. It is calculated from the exceedance probability $p$ — the probability that the event is exceeded in any single year:

T = \frac{1}{p} = \frac{1}{1 - F(x)}

(1)

where $F(x)$ is the cumulative distribution function (CDF) evaluated at intensity $x$ . For example, a 100-year return level has an exceedance probability of 1% per year.

We can invert the CDF to obtain the return level associated with a given return period: $x_T = F^{-1}(1 - 1/T)$ .

from plotly.subplots import make_subplots

# Distinct color per return period
rp_colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b']

# Compute return levels
return_periods = np.array([2, 5, 10, 20, 50, 100])
return_levels = genextreme.ppf(1 - 1/return_periods, *gev_params)

x = np.linspace(0, annual_max_df['rx1day'].max() + 25, 100)
pdf = genextreme.pdf(x, *gev_params)
cdf = genextreme.cdf(x, *gev_params)

# Build density histograms manually
counts_pdf, bin_edges = np.histogram(annual_max_df['rx1day'], bins=20, density=True)
bin_width = bin_edges[1] - bin_edges[0]
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
counts_cdf = np.cumsum(counts_pdf) * bin_width
y_max_pdf = max(counts_pdf.max(), pdf.max()) * 1.2

fig = make_subplots(rows=1, cols=2, subplot_titles=['PDF', 'CDF'])

# --- PDF panel (legend) ---
fig.add_trace(go.Bar(x=bin_centers, y=counts_pdf, width=bin_width * 0.95,
                     marker=dict(color='steelblue', line=dict(width=0)),
                     name='Annual maxima', legend='legend'), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=pdf, mode='lines', line=dict(color='indigo', width=2),
                         name='PDF of the fitted GEV', legend='legend'), row=1, col=1)

# --- CDF panel (legend2) ---
fig.add_trace(go.Bar(x=bin_centers, y=counts_cdf, width=bin_width * 0.95,
                     marker=dict(color='steelblue', line=dict(width=0)),
                     name='Annual maxima (cumulative)', legend='legend2'), row=1, col=2)
fig.add_trace(go.Scatter(x=x, y=cdf, mode='lines', line=dict(color='indigo', width=2),
                         name='CDF of the fitted GEV', legend='legend2'), row=1, col=2)

# --- Return level lines as traceable scatter traces (legend3) ---
for idx, (rp, rl) in enumerate(zip(return_periods, return_levels)):
    color = rp_colors[idx]
    hover_text = f'<b>{rp}-year return period</b><br>Return level: {rl:.1f} mm<br>Exceedance prob.: {1/rp:.1%}'

    # PDF panel — shown in legend3
    fig.add_trace(go.Scatter(
        x=[rl, rl], y=[0, y_max_pdf],
        mode='lines',
        line=dict(color=color, dash='dash', width=1.5),
        name=f'{rp}-yr ({1/rp:.0%})',
        legend='legend3', legendgroup=f'rp{rp}',
        showlegend=True,
        hovertemplate=hover_text + '<extra></extra>',
    ), row=1, col=1)

    # CDF panel — linked via legendgroup, hidden from legend
    fig.add_trace(go.Scatter(
        x=[rl, rl], y=[0, 1.05],
        mode='lines',
        line=dict(color=color, dash='dash', width=1.5),
        name=f'{rp}-yr ({1/rp:.0%})',
        legend='legend3', legendgroup=f'rp{rp}',
        showlegend=False,
        hovertemplate=hover_text + '<extra></extra>',
    ), row=1, col=2)

    # Horizontal annotation at top of each panel with background
    for col, xref, yref in [(1, 'x', 'y domain'), (2, 'x2', 'y2 domain')]:
        fig.add_annotation(
            x=rl, y=0.97, text=f'<b>{rp} yr</b>',
            showarrow=False,
            font=dict(size=9, color=color),
            xref=xref, yref=yref,
            bgcolor='rgba(255,255,255,0.85)',
            bordercolor=color, borderwidth=1, borderpad=2,
        )

fig.update_xaxes(title_text='Precipitation [mm]', showgrid=True, gridcolor='lightgrey')
fig.update_yaxes(title_text='Density', showgrid=True, gridcolor='lightgrey', row=1, col=1)
fig.update_yaxes(title_text='Cumulative Density', showgrid=True, gridcolor='lightgrey', row=1, col=2)
fig.update_layout(
    title=dict(
        text=f'Distributions of the Annual Maximum Daily Precipitation<br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].iloc[0]}, lat={point_df["lat"].iloc[0]}</sup>',
        font=dict(size=16),
    ),
    bargap=0,
    legend=dict(x=0.01, y=-0.15, xanchor='left', yanchor='top', font=dict(size=10),
                bgcolor='rgba(255,255,255,0.7)'),
    legend2=dict(x=0.51, y=-0.15, xanchor='left', yanchor='top', font=dict(size=10),
                 bgcolor='rgba(255,255,255,0.7)'),
    legend3=dict(x=0.01, y=-0.3, xanchor='left', yanchor='top', font=dict(size=10),
                 bgcolor='rgba(255,255,255,0.7)', orientation='h',
                 title=dict(text='Return periods:', font=dict(size=10))),
    width=1200,
    height=600,
)

fig.show()

Return level curves¶

To establish a clear relationship between the return period and the associated annual precipitation intensity (or return level), we can also plot these against each other.

return_periods = np.array([2, 5, 10, 20, 50, 100])
return_levels = genextreme.ppf(1 - 1/return_periods, *gev_params)

fig = go.Figure()

fig.add_trace(go.Scatter(
    x=return_periods,
    y=return_levels,
    mode='lines+markers+text',
    marker=dict(color='indigo', size=8),
    line=dict(color='indigo'),
    name='Return level curve from fitted GEV',
    text=[f'{rp} yr<br>{rl:.1f} mm' for rp, rl in zip(return_periods, return_levels)],
    textposition='top right',
    #mode='lines+markers+text',
))

fig.update_layout(
    #title=dict(text='Return level curve', font=dict(size=16)),
    title=dict(text=f'Return level curve<br><sup>{point_df['member_id'].unique()[0]}, lon={point_df['lon'].unique()[0]}, lat={point_df['lat'].unique()[0]}</sup>', font=dict(size=16)),
    xaxis=dict(title='Return period [years]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    yaxis=dict(title='Precipitation intensity [mm day⁻¹]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=0, y=1, xanchor='left', yanchor='top'),
    width=800,
    height=500,
)

fig.show()

# Print thresholds associated with return periods of 2, 5, 10, 20, 50, and 100 years
print("Precipitation intensity associated with different return periods at the selected point:")
for rp, rl in zip(return_periods, return_levels):
    print(f" - {rp}-year return period: {rl:.2f} mm year^-1")

Precipitation intensity associated with different return periods at the selected point:
 - 2-year return period: 38.87 mm year^-1
 - 5-year return period: 52.38 mm year^-1
 - 10-year return period: 61.14 mm year^-1
 - 20-year return period: 69.41 mm year^-1
 - 50-year return period: 79.92 mm year^-1
 - 100-year return period: 87.66 mm year^-1

Uncertainties embedded in the estimation of return periods¶

From the plots above, we can already define one uncertainty associated with the definition of return periods: the goodness of the fit.

However, there are also a few other uncertainties going into the estimation of the return period, that we briefly illustrate here.

Baseline period uncertainty¶

In the plots above, we used all data available for the historical model climate simulations.

print(f"The data used for fitting the GEV distribution covers {len(point_df['time'].dt.year.unique())} years and the period from {point_df['time'].min().date()} to {point_df['time'].max().date()}.")

The data used for fitting the GEV distribution covers 36 years and the period from 1970-01-01 to 2005-12-01.

The World Meteorological Organisation suggests using (at least) 30 years of data. So, let’s quickly assess the difference of using 30 years against the full time period.

# Extract data for baseline period defined in the settings above
print(f"\nCalculating return level curves for the baseline period ({baseline_start}-{baseline_end}) and comparing to the return level curve from the full data period.")
point_df_baseline = point_df[(point_df['year'] >= baseline_start) & (point_df['year'] <= baseline_end)]
gev_params_baseline = genextreme.fit(point_df_baseline['rx1day'])
return_levels_baseline = genextreme.ppf(1 - 1/return_periods, *gev_params_baseline)

print(f"{point_df['year'].min()}-{point_df['year'].max()}")
fig = go.Figure()

fig.add_trace(go.Scatter(
    x=return_periods,
    y=return_levels,
    mode='lines+markers+text',
    marker=dict(color='indigo', size=8),
    line=dict(color='indigo'),
    name=f'Full historical period ({point_df["year"].min()}-{point_df["year"].max()})',
    text=[f'{rp} yr<br>{rl:.1f} mm' for rp, rl in zip(return_periods, return_levels)],
    textposition='top right',
))

fig.add_trace(go.Scatter(
    x=return_periods,
    y=return_levels_baseline,
    mode='lines+markers+text',
    marker=dict(color='violet', size=8),
    line=dict(color='violet'),
    name=f'Baseline period ({baseline_start}-{baseline_end})',
    text=[f'{rp} yr<br>{rl:.1f} mm' for rp, rl in zip(return_periods, return_levels_baseline)],
    textposition='top right',
))

fig.update_layout(
    #title=dict(text='Return level curve', font=dict(size=16)),
    title=dict(text=f'Return level curve | Uncertainty associated with time period considered <br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].unique()[0]}, lat={point_df["lat"].unique()[0]}</sup>', font=dict(size=16)),
    xaxis=dict(title='Return period [years]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    yaxis=dict(title='Precipitation intensity [mm day⁻¹]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=0, y=1, xanchor='left', yanchor='top'),
    width=800,
    height=500,
)

fig.show()


Calculating return level curves for the baseline period (1976-2005) and comparing to the return level curve from the full data period.
1970-2005

Since fewer extremes are represented, the GEV based on a shorter time period is usually defining the return periods based on lower precipitation intensities.

Spatial uncertainty¶

Another uncertainty embedded in the estimation of the return periods is the spatial uncertainty. Climate models cannot resolve extremes at the exact location; so let’s use the extremes in a specific region (here NUTS2 region Molise, Italy) to illustrate how much the return level curve can vary across the region and how representative the single point is for the region as a whole.

# Use again the full data period but now plot the return level curve for all points in the region (i.e., for all lon/lat combinations) and compare to the return level curve from the single point used above. 
print(f"Read {sample_model_file}")
allpoints_df = pd.read_csv(sample_model_file, parse_dates=['year'])
allpoints_df['point_id'] = allpoints_df['lon'].astype(str) + '_' + allpoints_df['lat'].astype(str)
print(f"Number of unique points in the region: {allpoints_df['point_id'].nunique()}")
allpoints_df['year'] = allpoints_df['year'].dt.year

return_periods = np.array([2, 5, 10, 20, 50, 100])
return_levels_all_points = []
for point_id, group in allpoints_df.groupby('point_id'):
    annual_max_per_point = group.groupby('year')['rx1day'].max().dropna()
    if len(annual_max_per_point) < 5:
        continue
    gev_params = genextreme.fit(annual_max_per_point)
    return_levels_point = genextreme.ppf(1 - 1/return_periods, *gev_params)
    return_levels_all_points.append(return_levels_point)
return_levels_all_points = np.array(return_levels_all_points)
regional_mean_return_levels = return_levels_all_points.mean(axis=0)

# Plotly plot of the return level curves for all points in the region, with shading for the range of return levels across all points, and a line for the regional mean return level curve across all points, and a line for the return level curve from the single point used before above. Add hover information to show the return period and return level for each point, and add a legend to indicate which line corresponds to which point.
#import plotly.graph_objects as go   
fig = go.Figure()
# plot return level curves for all points in the region
for i in range(return_levels_all_points.shape[0]):  
    fig.add_trace(go.Scatter(x=return_periods, y=return_levels_all_points[i, :], mode='lines', line=dict(color='orchid', width=1), name='Single point GEV (all points in the region)' if i == 0 else None, showlegend=(i == 0)))
# add shading for the range of return levels across all points
fig.add_trace(go.Scatter(x=np.concatenate([return_periods, return_periods[::-1]]), y=np.concatenate([return_levels_all_points.min(axis=0), return_levels_all_points.max(axis=0)[::-1]]), fill='toself', fillcolor='rgba(255, 192, 203, 0.2)', line=dict(color='rgba(255, 192, 203, 0.2)'), 
                         name='Spatial uncertainty range (min–max)', showlegend=True))
# plot return level curve for the single point used before above
fig.add_trace(go.Scatter(x=return_periods, y=return_levels, mode='lines+markers', line=dict(color='indigo', width=2), marker=dict(size=8), name=f'Single point GEV (lon={point_df["lon"].unique()[0]}, lat={point_df["lat"].unique()[0]})', showlegend=True))
# plot regional mean return level curve across all points
fig.add_trace(go.Scatter(x=return_periods, y=regional_mean_return_levels, mode='lines+markers', line=dict(color='orchid', width=2), marker=dict(size=8), name='Mean return levels based on all points', showlegend=True))
fig.update_layout(title=f'Return level curve | Uncertainty associated with the spatial variability across the region<br><sup>{point_df["member_id"].unique()[0]}, region: {admin_id}</sup>', 
                  xaxis_title='Return period [years]', 
                  yaxis_title='Precipitation intensity [mm day<sup>-1</sup>]', 
                  legend=dict(x=0, y=1, xanchor='left', yanchor='top'),
                  width=800, height=500)
fig.show()

Read /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_amax_sample_model.csv
Number of unique points in the region: 45

print("Range of precipitation intensities across all points in the region for different return periods:")
for i, rp in enumerate(return_periods):
    print(f" - {rp}-year return period: {return_levels_all_points[:, i].min():.2f} to {return_levels_all_points[:, i].max():.2f} mm per day")

Range of precipitation intensities across all points in the region for different return periods:
 - 2-year return period: 34.66 to 57.66 mm per day
 - 5-year return period: 46.36 to 73.71 mm per day
 - 10-year return period: 51.02 to 84.96 mm per day
 - 20-year return period: 53.38 to 108.56 mm per day
 - 50-year return period: 55.46 to 150.64 mm per day
 - 100-year return period: 56.53 to 192.71 mm per day

Model uncertainty¶

So far, we have only used one climate model - but the power lies in using multiple models! So, let’s estimate the sensitivity of the return period to the climate model employed.

# read all models 
print(f"Reading {sample_point_file}")
sample_point_df = pd.read_csv(sample_point_file, parse_dates=['year'])
sample_point_df['year'] = sample_point_df['year'].dt.year
print(f"Number of unique ensemble members: {sample_point_df['member_id'].nunique()}")
unique_members = sample_point_df['member_id'].unique()

# use plotly to plot the return level curve for each member_id and remove the legend for each member_id to avoid overcrowding the legend, and add a single legend entry for all member_ids
import plotly.graph_objects as go
fig = go.Figure()
for member in unique_members:
    member_df = sample_point_df[sample_point_df['member_id'] == member]
    annual_max_per_member = member_df.groupby('year')['rx1day'].max().dropna()
    if len(annual_max_per_member) < 5:
        continue
    gev_params_member = genextreme.fit(annual_max_per_member)
    return_levels_member = genextreme.ppf(1 - 1/return_periods, *gev_params_member)
    fig.add_trace(go.Scatter(x=return_periods, y=return_levels_member, mode='lines+markers', name=f'{member}', line=dict(width=1), marker=dict(size=6)))
fig.update_layout(title=f'Return level curve | Uncertainty associated with the multi-model ensemble<br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].unique()[0]}, lat={point_df["lat"].unique()[0]}</sup>', 
                  xaxis_title='Return period [years]', 
                  yaxis_title='Precipitation intensity [mm day<sup>-1</sup>]', 
                  legend_title='Ensemble member', legend=dict(font=dict(size=10), bgcolor='rgba(255,255,255,0.5)'),
                  width=1200, height=600)
fig.show()

Reading /home/nejk/code/extreme_precip_exposure/data/ITF2/extreme_precip_exposure/ITF2_precip_hist_amax_sample_gridpoint.csv
Number of unique ensemble members: 71

Overview of uncertainty ranges¶

print("=" * 70)
print("UNCERTAINTY RANGES FOR THE CURRENT CASE")
print(f"Region: {admin_id} | Point: lon={point_df['lon'].iloc[0]}, lat={point_df['lat'].iloc[0]}")
print("=" * 70)

print(f"\n(i) Baseline period uncertainty  ({baseline_start}-{baseline_end} vs. full historical period)")
print(f"\n{'Return period':>16} | {'Full period':>12} | {'Baseline':>12} | {'Difference':>12}")
print("-" * 60)
for i, rp in enumerate(return_periods):
    full = return_levels[i]
    base = return_levels_baseline[i]
    print(f"{rp:>13}-year | {full:>10.1f} mm | {base:>10.1f} mm | {base-full:>+10.1f} mm")
print(f"\n  → Range of shift: {(return_levels_baseline - return_levels).min():+.1f} to {(return_levels_baseline - return_levels).max():+.1f} mm")

print(f"\n(ii) Spatial uncertainty  (range across all grid points in the region)")
print(f"\n{'Return period':>16} | {'Min':>9} | {'Max':>9} | {'Range':>9}")
print("-" * 52)
for i, rp in enumerate(return_periods):
    mn = return_levels_all_points[:, i].min()
    mx = return_levels_all_points[:, i].max()
    print(f"{rp:>13}-year | {mn:>7.1f} mm | {mx:>7.1f} mm | {mx-mn:>7.1f} mm")

print(f"\n(iii) Model uncertainty  (range across ensemble members at the selected point)")
_rl_members = []
for _member in sample_point_df['member_id'].unique():
    _mdf = sample_point_df[sample_point_df['member_id'] == _member]
    _amax = _mdf.groupby('year')['rx1day'].max().dropna()
    if len(_amax) >= 5:
        _p = genextreme.fit(_amax)
        _rl_members.append(genextreme.ppf(1 - 1/return_periods, *_p))
_rl_members = np.array(_rl_members)
print(f"\n{'Return period':>16} | {'Min':>9} | {'Max':>9} | {'Range':>9}")
print("-" * 52)
for i, rp in enumerate(return_periods):
    mn = _rl_members[:, i].min()
    mx = _rl_members[:, i].max()
    print(f"{rp:>13}-year | {mn:>7.1f} mm | {mx:>7.1f} mm | {mx-mn:>7.1f} mm")

print("\n" + "=" * 70)

======================================================================
UNCERTAINTY RANGES FOR THE CURRENT CASE
Region: ITF2 | Point: lon=14.5625, lat=41.6875
======================================================================

(i) Baseline period uncertainty  (1976-2005 vs. full historical period)

   Return period |  Full period |     Baseline |   Difference
------------------------------------------------------------
            2-year |       38.9 mm |       10.8 mm |      -28.1 mm
            5-year |       52.4 mm |       21.9 mm |      -30.5 mm
           10-year |       61.1 mm |       31.4 mm |      -29.8 mm
           20-year |       69.4 mm |       42.5 mm |      -26.9 mm
           50-year |       79.9 mm |       60.6 mm |      -19.3 mm
          100-year |       87.7 mm |       77.5 mm |      -10.1 mm

  → Range of shift: -30.5 to -10.1 mm

(ii) Spatial uncertainty  (range across all grid points in the region)

   Return period |       Min |       Max |     Range
----------------------------------------------------
            2-year |    34.7 mm |    57.7 mm |    23.0 mm
            5-year |    46.4 mm |    73.7 mm |    27.4 mm
           10-year |    51.0 mm |    85.0 mm |    33.9 mm
           20-year |    53.4 mm |   108.6 mm |    55.2 mm
           50-year |    55.5 mm |   150.6 mm |    95.2 mm
          100-year |    56.5 mm |   192.7 mm |   136.2 mm

(iii) Model uncertainty  (range across ensemble members at the selected point)

   Return period |       Min |       Max |     Range
----------------------------------------------------
            2-year |    24.1 mm |    64.2 mm |    40.1 mm
            5-year |    30.7 mm |    97.4 mm |    66.7 mm
           10-year |    34.9 mm |   127.4 mm |    92.5 mm
           20-year |    39.1 mm |   164.2 mm |   125.1 mm
           50-year |    44.6 mm |   226.9 mm |   182.3 mm
          100-year |    48.8 mm |   288.6 mm |   239.8 mm

======================================================================

fig = go.Figure()

# --- Spatial uncertainty band across all return periods ---
fig.add_trace(go.Scatter(
    x=return_periods, y=return_levels_all_points.min(axis=0),
    mode='lines', line=dict(color='orchid', width=0),
    showlegend=False, hoverinfo='skip',
))
fig.add_trace(go.Scatter(
    x=return_periods, y=return_levels_all_points.max(axis=0),
    mode='lines', line=dict(color='orchid', width=1.5),
    fill='tonexty', fillcolor='rgba(218, 112, 214, 0.2)',
    name='Spatial uncertainty (min–max)',
    hovertemplate='Spatial max: %{y:.1f} mm<extra></extra>',
))

# --- Model uncertainty band across all return periods ---
fig.add_trace(go.Scatter(
    x=return_periods, y=_rl_members.min(axis=0),
    mode='lines', line=dict(color='lightseagreen', width=0),
    showlegend=False, hoverinfo='skip',
))
fig.add_trace(go.Scatter(
    x=return_periods, y=_rl_members.max(axis=0),
    mode='lines', line=dict(color='lightseagreen', width=1.5),
    fill='tonexty', fillcolor='rgba(32, 178, 170, 0.2)',
    name='Model uncertainty (min–max)',
    hovertemplate='Model max: %{y:.1f} mm<extra></extra>',
))

# --- Full historical period return level curve ---
fig.add_trace(go.Scatter(
    x=return_periods, y=return_levels,
    mode='lines+markers', line=dict(color='indigo', width=2.5),
    marker=dict(size=7),
    name='Full historical period',
    hovertemplate='%{x}-yr return level (full period): %{y:.1f} mm<extra></extra>',
))

# --- Baseline period return level curve ---
fig.add_trace(go.Scatter(
    x=return_periods, y=return_levels_baseline,
    mode='lines+markers', line=dict(color='mediumpurple', width=2.5, dash='dash'),
    marker=dict(size=7),
    name=f'Baseline period ({baseline_start}–{baseline_end})',
    hovertemplate='%{x}-yr return level (baseline): %{y:.1f} mm<extra></extra>',
))

fig.update_layout(
    title=f'Return level curves and uncertainty ranges<br><sup>{point_df["member_id"].unique()[0]}, lon={point_df["lon"].unique()[0]}, lat={point_df["lat"].unique()[0]}</sup>',
    xaxis=dict(
        title='Return period [years]',
        tickmode='array',
        tickvals=return_periods.tolist(),
        ticktext=[f'{rp} yr' for rp in return_periods],
        showgrid=True, gridcolor='lightgrey', gridwidth=0.5,
    ),
    yaxis=dict(title='Precipitation intensity [mm day⁻¹]', showgrid=True, gridcolor='lightgrey', gridwidth=0.5),
    legend=dict(x=0, y=1, xanchor='left', yanchor='top'),
    width=900,
    height=500,
)
fig.show()

Summary¶

How a return period is estimated¶

The return period $T$ of an extreme precipitation event is defined as the average time between two events of the same or greater intensity. It is calculated from the exceedance probability

In practice, the estimation proceeds in three steps:

Calculate annual maxima from the, e.g., daily precipitation timeseries.
Fit a Generalised Extreme Value (GEV) distribution to those annual maxima.
Invert the CDF to obtain the return level associated with a given return period: $x_T = F^{-1}(1 - 1/T)$ .

The return level curve plots $x_T$ against $T$ , providing a compact view of how extreme precipitation intensities scale with rarity.

Sources of uncertainty¶

The return period estimate carries several layers of uncertainty:

(i) Goodness of fit of the GEV distribution
The GEV is a theoretical model. If the annual maxima do not follow a GEV well — e.g. due to non-stationarity, seasonality, or mixed precipitation regimes, and the fitted parameters and resulting return levels will be biased. Diagnostics such as Q-Q plots and the PDF/CDF overlays shown above help assess this.

(ii) Baseline period length
The length of the record used to fit the GEV directly affects the stability of the estimated parameters. Shorter records mean fewer annual maxima, leading to larger estimation uncertainty, particularly for high return periods (e.g. 50- or 100-year events) that are poorly sampled in short records. The WMO recommends a minimum of 30 years.

(iii) Spatial uncertainty
Climate model output is provided on a spatial grid. The return level at any single grid point may not represent the broader region. Across the grid points in the study region, return levels can vary substantially, reflecting real spatial gradients in precipitation extremes as well as model artefacts.

(iv) Model uncertainty
Different climate models, driven by different parameterisations and physical assumptions, can produce substantially different precipitation extremes, even when driven by the same boundary conditions. The multi-model ensemble spread provides a measure of this structural uncertainty.

Take-aways and best practice

A single grid point from a single climate model is not sufficient to robustly estimate return periods. The uncertainty ranges demonstrated above show that both the choice of model and the choice of location can shift return level estimates by tens of millimetres — often by more than the signal itself at shorter return periods.

Best practice recommendations:

Use multiple models. Always derive return levels from a climate model ensemble rather than a single model run. Report the full spread as a measure of model uncertainty.
Account for spatial variability. Return levels vary across the region. Rather than relying on a single grid point, compute return levels at all grid points within your area of interest and consider the regional distribution (mean, range, percentiles).
Use an adequate record length. At least 30 years of annual maxima are needed to fit a GEV reliably. Shorter records systematically underestimate return levels for rare events.
Be transparent about uncertainty. When communicating a return level (e.g. “the 100-year daily precipitation is X mm”), always accompany it with the uncertainty range across models and grid points.