docs/notebooks/11_Climatological_ESP_forecasting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11 - Climatological ESP forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extended Streamflow Prediction (ESP) forecasts from climatological time series\n",
    "\n",
    "This notebook shows how to perform a climatological Extended Streamflow Prediction (ESP) forecast, using historical weather as a proxy for future weather.\n",
    "\n",
    "The general idea is to initialize the state of the hydrological model to represent current conditions, but instead of using weather forecasts to predict future flows, we run the model with observed, historical weather series from past years. So for example if we have 30 years of weather observations, we get 30 different forecasts. The accuracy of this forecast ensemble can then be evaluated by different probabilistic metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from ravenpy.config import commands as rc\n",
    "from ravenpy.config.emulators import GR4JCN\n",
    "from ravenpy.utilities import forecasting\n",
    "from ravenpy.utilities.testdata import get_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the model simulations\n",
    "\n",
    "Here we set model parameters somewhat arbitrarily, but you can set the parameters to the calibrated parameters as seen in the \"06_Raven_calibration\" notebook we previously encountered.\n",
    "\n",
    "We also need to choose the forecast issue date. Each forecast will start with the same day and month. For example, jun-06-1980 will compare the climatology using all jun-06's from the dataset. Finally, we can provide the forecast duration (in number of days) as well as the historical meteorological years we want to use to generate the ESP forecast. This allows selecting years that we want to include in the forecast. For example, perhaps we only want to generate a forecast using wet or dry years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the selected watershed's time series. You can use your own time-series for your catchment by replacing\n",
    "# this line with the name / path of your input file.\n",
    "ts = get_file(\"raven-gr4j-cemaneige/Salmon-River-Near-Prince-George_meteo_daily.nc\")\n",
    "\n",
    "# This is the forecast start date, on which the forecasts will be launched.\n",
    "start_date = dt.datetime(1980, 6, 1)\n",
    "\n",
    "# Provide the length of the forecast, in days:\n",
    "forecast_duration = 100\n",
    "\n",
    "# Define HRU to build the hydrological model\n",
    "hru = {}\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",
    "            ts, data_type=data_type, alt_names=alt_names, data_kwds=data_kwds\n",
    "        )\n",
    "    ],\n",
    "    HRUs=[hru],\n",
    "    StartDate=start_date,\n",
    "    Duration=forecast_duration,\n",
    "    RunName=\"full\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Issuing the ESP forecast\n",
    "\n",
    "Here we launch the code that will perform the ESP forecast. Depending on the number of years in the historical dataset and the forecast duration, it might take a while to return a forecast result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate the climatological ESP:\n",
    "ESP_sims = forecasting.climatology_esp(\n",
    "    config=model_config,\n",
    "    years=[\n",
    "        1982,\n",
    "        1998,\n",
    "        2003,\n",
    "        2004,\n",
    "    ],  # List of years to use in the forecast. Optional. Will use all years by default.\n",
    ")\n",
    "\n",
    "# Show the results in an xarray dataset, ready to use:\n",
    "ESP_sims.hydrograph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now inspect and graph the resulting climatological ESP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ESP_sims.hydrograph.q_sim[:, :, 0].plot.line(x=\"time\")\n",
    "plt.title(\"GR4JCN climatological ESP for 1980-06-01\")\n",
    "plt.xticks(rotation=90)\n",
    "plt.grid(\"on\")\n",
    "plt.xlabel(\"Time [days]\")\n",
    "plt.ylabel(\"Streamflow $[m^3s^{-1}]$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the forecast scores\n",
    "\n",
    "There are different metric to evaluate the performance of forecasts. As an example, here we are computing the CRPS metric, using the [xskillscore](https://xskillscore.readthedocs.io/en/stable/) library included in PAVICS-Hydro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import xskillscore as xs\n",
    "\n",
    "# Align time axes to get the observed streamflow time series for the same time frame as the ESP forecast ensemble\n",
    "q_obs, q_sims = xr.align(xr.open_dataset(ts).qobs, ESP_sims.hydrograph, join=\"inner\")\n",
    "\n",
    "# Adjust the streamflow to convert missing data from -1.2345 format to NaN. Set all negative values to NaN.\n",
    "q_obs = q_obs.where(q_obs > 0, np.nan)\n",
    "\n",
    "# Compute the Continuous Ranked Probability Score using xskillscore\n",
    "xs.crps_ensemble(q_obs, q_sims, dim=\"time\").q_sim.values[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing a climatology ESP hindcast\n",
    "In this section, we make the hindcasts for each initialization date that we desire. Here we will extract ESP forecasts for a given calendar date for the years in \"hindcast_years\" as hindcast dates. Each ESP hindcast uses all available data in the `ts` dataset, so in this case we will have 56/57 members for each hindcast initialization depending on the date that we start on, UNLESS we specify a list of years manually. The \"hindcasts\" dataset generated contains all of the flow data from the ESP hindcasts for the initialization dates. The `q_obs` dataset contains all q_obs in the timeseries: Climpred will sort it all out during its processing. Note that the format of these datasets is tailor-made to be used in climpred, and thus has specific dimension names.\n",
    "\n",
    "This is a slimmed down example of how we would run an ESP forecast over multiple years to assess the skill of such a forecast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hindcasts = forecasting.hindcast_climatology_esp(\n",
    "    config=model_config,  # Note that the forecast duration is already set-up in the model_config above.\n",
    "    warm_up_duration=365,  # number of days for the warm-up\n",
    "    years=[1985, 1986, 1987, 1988, 1989, 1990],\n",
    "    hindcast_years=[2001, 2002, 2003, 2004, 2005, 2006, 2007],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the forecast using different metrics\n",
    "Once we have the correctly formatted datasets, Make the hindcast object for climpred\n",
    "\n",
    "These three functions respectively compute the rank histogram, the CRPS and the reliability for the set of initialized dates (i.e. forecast issue dates, here 1 day per year at the same calendar day)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once we have the correctly formatted datasets, Make the hindcast object for climpred\n",
    "\n",
    "# We first need to get the observed streamflow:\n",
    "q_obs = xr.open_dataset(ts)\n",
    "\n",
    "# However, our simulated streamflow is named \"q_sim\" and climpred requires the observation to be named the same thing\n",
    "# so let's rename it. While we're at it, we need to make sure that the identifier is the same. In our observation\n",
    "# dataset, it is called \"nstations\" but in our simulated streamflow it's called \"nbasins\". Here we standardize.\n",
    "q_obs = q_obs.rename({\"qobs\": \"q_sim\", \"nstations\": \"nbasins\"})\n",
    "\n",
    "# Make the hindcasting object we can use to compute statistics and metrics\n",
    "hindcast_object = forecasting.to_climpred_hindcast_ensemble(hindcasts, q_obs)\n",
    "\n",
    "\n",
    "# This function is used to convert to binary to see if yes/no forecast is larger than observations\n",
    "def pos(x):\n",
    "    return x > 0  # Check for binary outcome\n",
    "\n",
    "\n",
    "# Rank histogram verification metric\n",
    "rank_histo_verif = hindcast_object.verify(\n",
    "    metric=\"rank_histogram\",\n",
    "    comparison=\"m2o\",\n",
    "    dim=[\"member\", \"init\"],\n",
    "    alignment=\"same_inits\",\n",
    ")\n",
    "# CRPS verification metric\n",
    "crps_verif = hindcast_object.verify(\n",
    "    metric=\"crps\",\n",
    "    comparison=\"m2o\",\n",
    "    dim=[\"member\", \"init\"],\n",
    "    alignment=\"same_inits\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can explore and plot the CRPS as a function of lead-time, for example. Results are stored as a dataset and\n",
    "# can thus be integrated into any simulation or processes.\n",
    "plt.plot(crps_verif.q_sim)\n",
    "plt.xlabel(\"Lead time [days]\")\n",
    "plt.ylabel(\"CRPS $[m^3s^{-1}]$\")\n",
    "plt.grid(\"on\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}