Probabilistic flood risk assessment

In this notebook, we combine the forecasting abilities and the time series analysis capabilities in a single seamless process to estimate the flood risk of a probabilistic forecast. As an example, we first perform a frequency analysis on an observed time series, then estimate the streamflow associated to a 2-year return period. We then perform a climatological ESP forecast (to ensure repeatability, but a realtime forecast would work too!) and estimate the probability of flooding (exceeding the threshold) given the ensemble of members in the probabilistic forecast.

import warnings

from numba.core.errors import NumbaDeprecationWarning

warnings.simplefilter("ignore", category=NumbaDeprecationWarning)

%matplotlib inline

import datetime as dt

import xclim
from matplotlib import pyplot as plt

from ravenpy.utilities.testdata import get_file, open_dataset

Perform the time series analysis on observed data for the catchment using the frequency analysis WPS capabilities.

# Get the data that we will be using for the demonstration.
file = "raven-gr4j-cemaneige/Salmon-River-Near-Prince-George_meteo_daily.nc"
ts = open_dataset(file).qobs

# Perform the frequency analysis for various return periods. We compute 2, 5, 10, 25, 50 and 100 year return
# periods, but later on we will only compare the forecasts to the 2 year return period.
out = xclim.generic.return_level(
    ts, mode="max", t=(2, 5, 10, 25, 50, 100), dist="gumbel_r"
)
out

# Plot the results of the flows as a function of return period.
fig, ax = plt.subplots(1)
lines = out.plot(ax=ax)

# Get 2-year return period from the frequency analysis
threshold = out.sel(return_period=2).values
print(f"Threshold: {threshold:.1f}")

pt = ax.plot([2], [threshold], "ro")

ax.annotate(
    "Flow threshold, set at 2-year return period",
    (2, threshold),
    xytext=(25, 10),
    textcoords="offset points",
    arrowprops=dict(arrowstyle="->", connectionstyle="arc3"),
)

Probabilistic forecast

In this example, we will perform an ensemble hydrological forecast and will then compute the probability of flooding given a flooding threshold. Start by building the model configuration as in the Tutorial Notebook 11:

from ravenpy.config import commands as rc
from ravenpy.config.emulators import GR4JCN
from ravenpy.utilities.forecasting import climatology_esp, compute_forecast_flood_risk

# Choose the forecast date. Each forecast will start with the same day and month.
# For example, jan-05-2001 will compare the climatology using all jan-05ths from the dataset)
fdate = dt.datetime(2003, 4, 13)

# The dataset to use to get the forecast timeseries:
duration = 30  # Length in days of the climatological ESP forecast

# Define HRU to build the hydrological model
hru = dict(
    area=4250.6,
    elevation=843.0,
    latitude=54.4848,
    longitude=-123.3659,
    hru_type="land",
)

# Set alternative names for netCDF variables
alt_names = {
    "TEMP_MIN": "tmin",
    "TEMP_MAX": "tmax",
    "RAINFALL": "rain",
    "SNOWFALL": "snow",
}

# Data types to extract from netCDF
data_type = ["TEMP_MAX", "TEMP_MIN", "RAINFALL", "SNOWFALL"]
data_kwds = {
    "ALL": {
        "elevation": hru[
            "elevation"
        ],  # No need for lat/lon as they are included in the netcdf file already
    }
}
# Model configuration
model_config = GR4JCN(
    params=[0.529, -3.396, 407.29, 1.072, 16.9, 0.947],
    Gauge=[
        rc.Gauge.from_nc(
            get_file(file),
            data_type=data_type,
            alt_names=alt_names,
            data_kwds=data_kwds,
        )
    ],
    HRUs=[hru],
    StartDate=fdate,
    Duration=duration,
    RunName="Probabilistic_flood_risk_NB",
)

Now that the configuration is ready, launch the ESP forecasting tool to generate an ensemble hydrological forecast:

# Launch the ESP forecasting method
ESP_sims = climatology_esp(
    config=model_config,
)

# Show the results in an xarray dataset, ready to use:
ESP_sims.hydrograph

# Plot the forecasts and the 2-year threshold previously estimated
fig, ax = plt.subplots(1)
ESP_sims.hydrograph.q_sim[:, :, 0].plot.line(
    ax=ax, hue="member", add_legend=False, color="gray", lw=0.5
)
t = ax.axhline(threshold, color="red")

# Now compute the flood risk given the probabilistic forecast and the threshold associated to the 2-year return
# period.

threshold = out.sel(return_period=2).values

# Run the flood forecast risk tool to extract the probability of exceedance in netcdf format and xarray Dataset format
flood_risk_data = compute_forecast_flood_risk(
    forecast=ESP_sims.hydrograph.q_sim,
    flood_level=threshold,
)

# Extract the data and plot
fig, ax = plt.subplots(1)
l = flood_risk_data.exceedance_probability.plot()
ax.set_ylabel("Flood risk")

Results analysis

We can see from the above figure that there is no risk of exceeding the 2-year return period for the selected dates of the forecast.