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.

%matplotlib inline
import datetime as dt
import warnings

import xarray as xr
import xclim
from matplotlib import pyplot as plt
from numba.core.errors import NumbaDeprecationWarning

# Utility that simplifies working with test data hosted on GitHub
from ravenpy.testing.utils import get_file

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

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 = xr.open_dataset(get_file(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 import emulators
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 = emulators.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.