docs/notebooks/Assess_probabilistic_flood_risk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Probabilistic flood risk assessment\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "from numba.core.errors import NumbaDeprecationWarning\n",
    "\n",
    "warnings.simplefilter(\"ignore\", category=NumbaDeprecationWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import datetime as dt\n",
    "\n",
    "import xclim\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from ravenpy.utilities.testdata import get_file, open_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "Perform the time series analysis on observed data for the catchment using the frequency analysis WPS capabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the data that we will be using for the demonstration.\n",
    "file = \"raven-gr4j-cemaneige/Salmon-River-Near-Prince-George_meteo_daily.nc\"\n",
    "ts = open_dataset(file).qobs\n",
    "\n",
    "# Perform the frequency analysis for various return periods. We compute 2, 5, 10, 25, 50 and 100 year return\n",
    "# periods, but later on we will only compare the forecasts to the 2 year return period.\n",
    "out = xclim.generic.return_level(\n",
    "    ts, mode=\"max\", t=(2, 5, 10, 25, 50, 100), dist=\"gumbel_r\"\n",
    ")\n",
    "out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results of the flows as a function of return period.\n",
    "fig, ax = plt.subplots(1)\n",
    "lines = out.plot(ax=ax)\n",
    "\n",
    "# Get 2-year return period from the frequency analysis\n",
    "threshold = out.sel(return_period=2).values\n",
    "print(f\"Threshold: {threshold:.1f}\")\n",
    "\n",
    "pt = ax.plot([2], [threshold], \"ro\")\n",
    "\n",
    "ax.annotate(\n",
    "    \"Flow threshold, set at 2-year return period\",\n",
    "    (2, threshold),\n",
    "    xytext=(25, 10),\n",
    "    textcoords=\"offset points\",\n",
    "    arrowprops=dict(arrowstyle=\"->\", connectionstyle=\"arc3\"),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## Probabilistic forecast\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "from ravenpy.config import commands as rc\n",
    "from ravenpy.config.emulators import GR4JCN\n",
    "from ravenpy.utilities.forecasting import climatology_esp, compute_forecast_flood_risk\n",
    "\n",
    "# Choose the forecast date. Each forecast will start with the same day and month.\n",
    "# For example, jan-05-2001 will compare the climatology using all jan-05ths from the dataset)\n",
    "fdate = dt.datetime(2003, 4, 13)\n",
    "\n",
    "# The dataset to use to get the forecast timeseries:\n",
    "duration = 30  # Length in days of the climatological ESP forecast\n",
    "\n",
    "# Define HRU to build the hydrological model\n",
    "hru = dict(\n",
    "    area=4250.6,\n",
    "    elevation=843.0,\n",
    "    latitude=54.4848,\n",
    "    longitude=-123.3659,\n",
    "    hru_type=\"land\",\n",
    ")\n",
    "\n",
    "# Set alternative names for netCDF variables\n",
    "alt_names = {\n",
    "    \"TEMP_MIN\": \"tmin\",\n",
    "    \"TEMP_MAX\": \"tmax\",\n",
    "    \"RAINFALL\": \"rain\",\n",
    "    \"SNOWFALL\": \"snow\",\n",
    "}\n",
    "\n",
    "# Data types to extract from netCDF\n",
    "data_type = [\"TEMP_MAX\", \"TEMP_MIN\", \"RAINFALL\", \"SNOWFALL\"]\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\n",
    "            \"elevation\"\n",
    "        ],  # No need for lat/lon as they are included in the netcdf file already\n",
    "    }\n",
    "}\n",
    "# Model configuration\n",
    "model_config = GR4JCN(\n",
    "    params=[0.529, -3.396, 407.29, 1.072, 16.9, 0.947],\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            get_file(file),\n",
    "            data_type=data_type,\n",
    "            alt_names=alt_names,\n",
    "            data_kwds=data_kwds,\n",
    "        )\n",
    "    ],\n",
    "    HRUs=[hru],\n",
    "    StartDate=fdate,\n",
    "    Duration=duration,\n",
    "    RunName=\"Probabilistic_flood_risk_NB\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "Now that the configuration is ready, launch the ESP forecasting tool to generate an ensemble hydrological forecast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Launch the ESP forecasting method\n",
    "ESP_sims = climatology_esp(\n",
    "    config=model_config,\n",
    ")\n",
    "\n",
    "# Show the results in an xarray dataset, ready to use:\n",
    "ESP_sims.hydrograph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Plot the forecasts and the 2-year threshold previously estimated\n",
    "fig, ax = plt.subplots(1)\n",
    "ESP_sims.hydrograph.q_sim[:, :, 0].plot.line(\n",
    "    ax=ax, hue=\"member\", add_legend=False, color=\"gray\", lw=0.5\n",
    ")\n",
    "t = ax.axhline(threshold, color=\"red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Now compute the flood risk given the probabilistic forecast and the threshold associated to the 2-year return\n",
    "# period.\n",
    "\n",
    "threshold = out.sel(return_period=2).values\n",
    "\n",
    "# Run the flood forecast risk tool to extract the probability of exceedance in netcdf format and xarray Dataset format\n",
    "flood_risk_data = compute_forecast_flood_risk(\n",
    "    forecast=ESP_sims.hydrograph.q_sim,\n",
    "    flood_level=threshold,\n",
    ")\n",
    "\n",
    "# Extract the data and plot\n",
    "fig, ax = plt.subplots(1)\n",
    "l = flood_risk_data.exceedance_probability.plot()\n",
    "ax.set_ylabel(\"Flood risk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "### Results analysis\n",
    "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.\n"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}