Joint Probability Example

Joint Probability Example#

This notebook demonstrates how to compute the probability of joint water level and precipitation extremes.
The analysis considers two cases:
- wl_prcond: Water level peaks-over-threshold events paired with the maximum precipitation event occurring within ±1 day.
- pr_wlcond: Precipitation peaks-over-threshold events paired with the maximum water level event occurring within ±1 day.
Univariate distributions are fitted to the marginals (water levels and precipitation) separately. Joint events are then analyzed using copulas, which are selected based on the Bayesian Information Criterion (BIC).
For methodology and limitations, see the technical report: Braschi, L.C., & Leys, V. (2025). Extrêmes de niveaux d’eau côtiers et extrêmes conjoints de niveaux d’eau et de précipitations au Canada. Rapport présenté à Ouranos. Montréal: CBCL Limited, 37 p.

The data for this notebook is stored in the cloud, so the first step is to access it.

Imports and Paths#

import pandas as pd
from scipy.stats import kendalltau
import xarray as xr
from pathlib import Path
import s3fs


from peach.frontend.wl_parameters import IndicatorObsWL, IndicatorSimWL
from peach.frontend.cbcl_utils import matching_events
from peach.frontend.jp_parameters import (
    IndicatorObsWLCOND,
    IndicatorObsPRPOT,
    IndicatorObsPRCOND,
    IndicatorSimPRPOT,
    IndicatorSimWLCOND,
    IndicatorSimPRPOT,
    IndicatorSimPRCOND,
    IndicatorObsJP,
    IndicatorSimJP
)

copulas = ["gaussian", "student", "clayton", "frank", "gumbel"]
scipy_dists = ["norm", "t", "gamma", 'genextreme']

stn_num = '07330'
ahccd_num = '1018620'

s3r = s3fs.S3FileSystem(anon=True, use_ssl=False, client_kwargs={"endpoint_url": "https://minio.ouranos.ca"})

def open_link(filename: str) -> str:
    """Read NetCDF object from MinIO server and return as a string."""
    root = "portail-ing/tutorial_data"
    link = s3r.open(f"{root}/{filename}")
    return link

# Peaks-over-threshold (pot) data
wl_pot_file = f"{stn_num}_wl_pot.nc"
pr_pot_file = f"{ahccd_num}_pr_pot.nc"

# Timeseries data (used to select conditional extremes)
wl_timeseries_file = f"{stn_num}_wl.nc"
pr_timeseries_file = f"{ahccd_num}_pr.nc"

# Future simulations: decadal sea-level delta and daily precipitation timeseries
sl_sim_file = f"{stn_num}_sl.nc"
pr_sim_file = f"{ahccd_num}_pr_sim.nc"

WARNING:matplotlib.font_manager:Matplotlib is building the font cache; this may take a moment.

config.py - Using local cache: True, with file None

2026-02-09 20:10:53,324 | INFO | peach.frontend.parameters | Logger successfully configured

2026-02-09 20:10:53,366 | INFO | peach.frontend.parameters | Workspace: /home/docs/checkouts/readthedocs.org/user_builds/cs-peach/conda/latest/lib/python3.14/workspace

2026-02-09 20:10:53,368 | INFO | URLParameterized | Logger successfully configured

2026-02-09 20:10:53,381 | INFO | peach.frontend.idf_parameters | Logger successfully configured

Step 1 - User Settings#

ref_period = (1995, 2014)
fut_period = (2070, 2100)
pval_thresh = 0.05 # p-value for the kendall dependence test
pr_rp = 10 # Return period for precipitation
wl_rp = 10 # Return period for water level

Step 2 - Load Data#

# Load peaks-over-threshold (pot) data
wl_pot_backend = xr.open_dataarray(open_link(wl_pot_file))
pr_pot_backend = xr.open_dataarray(open_link(pr_pot_file))

# Load timeseries data
wl = xr.open_dataarray(open_link(wl_timeseries_file))
pr = xr.open_dataarray(open_link(pr_timeseries_file))

# Load and format future simulation data
sl_sim_backend = xr.open_dataset(open_link(sl_sim_file))['sl_delta']
sim = xr.open_dataarray(open_link(pr_sim_file))
multi_index = pd.MultiIndex.from_arrays(
    [sim['variant_label'].values, sim['source_id'].values, sim['experiment_id'].values],
    names=["variant_label", "source_id", "experiment_id"]
)
pr_sim_backend = sim.assign_coords(realization=("realization", multi_index))

Step 3 - Conditional Events#

# Select maximum precipitation events occuring within 1 day of water level peaks
wl_pot_backend, pr_cond_backend = matching_events(pot_da=wl_pot_backend, timeseries_da=pr)

# Select maximum water levels events occuring within 1 day of precpitation peaks
pr_pot_backend, wl_cond_backend = matching_events(pot_da=pr_pot_backend, timeseries_da=wl)

Step 4 - Dependence#

# Assess statistical dependence with Kendall's Tau
wl_prcond_ktau, wl_prcond_pval = kendalltau(wl_pot_backend.values, pr_cond_backend.values)
pr_wlcond_ktau, pr_wlcond_pval = kendalltau(pr_pot_backend.values, wl_cond_backend.values)

def sig_string(pval, pval_thresh):
    if pval <= pval_thresh:
        return 'SIGNIFICANT'
    else:
        return 'NOT SIGNIFICANT'

print(f'Water level POT & conditional precipitation extremes: {sig_string(wl_prcond_pval, pval_thresh)}\n'
      f'p-value {round(wl_prcond_pval, 4)} and kendall tau {round(wl_prcond_ktau, 2)}')
print(f'Precipitation POT & conditional water level extremes: {sig_string(pr_wlcond_pval, pval_thresh)}\n'
      f'p-value {round(pr_wlcond_pval, 4)} and kendall tau {round(pr_wlcond_ktau, 2)}')

Water level POT & conditional precipitation extremes: SIGNIFICANT
p-value 0.0 and kendall tau 0.18
Precipitation POT & conditional water level extremes: SIGNIFICANT
p-value 0.0 and kendall tau 0.15

Step 5 - Joint analysis for wl_prcond#

# Marginals for observations
wl_pot = IndicatorObsWL(data=wl_pot_backend, period=ref_period)
pr_cond = IndicatorObsPRCOND(data=pr_cond_backend, period=ref_period)

# Marginals for simulations
wl_sim = IndicatorSimWL(obs=wl_pot, data=sl_sim_backend, period=fut_period)
pr_sim = IndicatorSimPRCOND(obs=pr_cond, data=pr_sim_backend, wl_pot=wl_pot_backend, period=fut_period)

# Copula for observations
jp_obs = IndicatorObsJP(obs_pot=wl_pot, obs_cond=pr_cond, name="wl_prcond", period=ref_period)

# Copula for simulations
jp_sim = IndicatorSimJP(sim_pot=wl_sim, sim_cond=pr_sim, obs_cop=jp_obs, period=fut_period)

# Example joint exceedance assessment for wl_prcond
result_obs = round(jp_obs.sf([wl_pot.isf(1/wl_rp).item(), pr_cond.isf(1/pr_rp).item()]).item(),2)
result_sim = round(jp_sim.sf([wl_pot.isf(1/wl_rp).item(), pr_cond.isf(1/pr_rp).item()]).item(),2)

result_text = (f"""
    The probability of precipitation exceeding a 1-in-{pr_rp} event and water levels exceeding a 1-in-{wl_rp} event
    during the reference period {ref_period} is {round((1/pr_rp)*(1/wl_rp), 4)} under the assumption of independence
    (0.1 x 0.1 = 0.01).

    Under the present analysis, the estimated probability of joint exceedance is {result_obs} for the reference period
    {ref_period} and {result_sim} for the future period {fut_period}.
    """)
print(result_text)

    The probability of precipitation exceeding a 1-in-10 event and water levels exceeding a 1-in-10 event
    during the reference period (1995, 2014) is 0.01 under the assumption of independence
    (0.1 x 0.1 = 0.01).

    Under the present analysis, the estimated probability of joint exceedance is 0.03 for the reference period
    (1995, 2014) and 0.43 for the future period (2070, 2100).
    

Step 6 - Joint analysis for pr_wlcond#

# Marginals for observations
pr_pot = IndicatorObsPRPOT(data=pr_pot_backend, period=ref_period)
wl_cond = IndicatorObsWLCOND(data=wl_cond_backend, period=ref_period)

# Marginals for simulations
pr_sim = IndicatorSimPRPOT(obs=pr_pot, data=pr_sim_backend, period=fut_period)
wl_sim = IndicatorSimWLCOND(obs=wl_cond, data=sl_sim_backend, period=fut_period)

# Copula for observations
jp_obs = IndicatorObsJP(obs_pot=pr_pot, obs_cond=wl_cond, name="pr_wlcond", period=ref_period)

# Copula for simulations
jp_sim = IndicatorSimJP(sim_pot=pr_sim, sim_cond=wl_sim, obs_cop=jp_obs, period=fut_period)

# Example joint exceedance assessment for pr_wlcond
result_obs = round(jp_obs.sf([pr_pot.isf(1/pr_rp).item(), wl_cond.isf(1/wl_rp).item()]).item(),2)
result_sim = round(jp_sim.sf([pr_pot.isf(1/pr_rp).item(), wl_cond.isf(1/wl_rp).item()]).item(),2)

result_text = (f"""
    The probability of precipitation exceeding a 1-in-{pr_rp} event and water levels exceeding a 1-in-{wl_rp} event
    during the reference period {ref_period} is {round((1/pr_rp)*(1/wl_rp), 4)} under the assumption of independence
    (0.1 x 0.1 = 0.01).

    Under the present analysis, the estimated probability of joint exceedance is {result_obs} for the reference period
    {ref_period} and {result_sim} for the future period {fut_period}.
    """)
print(result_text)

    The probability of precipitation exceeding a 1-in-10 event and water levels exceeding a 1-in-10 event
    during the reference period (1995, 2014) is 0.01 under the assumption of independence
    (0.1 x 0.1 = 0.01).

    Under the present analysis, the estimated probability of joint exceedance is 0.05 for the reference period
    (1995, 2014) and 0.12 for the future period (2070, 2100).
    

# The reliability of these results should be assessed based on confidence in event selection, marginal fits, and copula fit.

# Use the `dist` method to see which copula was selected and the `bic` method to see its Bayesian Information Criterion (BIC).
print(f'The BIC for the {jp_obs.dist} copula is: {round(jp_obs.bic.item(),4)}')

The BIC for the gumbel copula is: -0.0517

Joint Probability Example

Contents

Joint Probability Example#

Imports and Paths#

Step 1 - User Settings#

Step 2 - Load Data#

Step 3 - Conditional Events#

Step 4 - Dependence#

Step 5 - Joint analysis for wl_prcond#

Step 6 - Joint analysis for pr_wlcond#