docs/notebooks/Comparing_hindcasts_and_ESP_forecasts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing hindcasts to a climatological ensemble streamflow prediction (ESP)\n",
    "\n",
    "This notebook shows how to use climatological weather to perform a Climatology-based Extended Streamflow Prediction (ESP) forecast. Then using the same initial states, uses the CaSPar archived weather forecasts to generate streamflow hindcasts over the same period. It is thus possible to compare both approaches.\n",
    "\n",
    "CaSPAr (Canadian Surface Prediction Archive) is an archive of historical ECCC forecasts developed by Juliane Mai at the University of Waterloo, Canada. More details on CaSPAr can be found here https://caspar-data.ca/.\n",
    "\n",
    "Mai, J., Kornelsen, K.C., Tolson, B.A., Fortin, V., Gasset, N., Bouhemhem, D., Schäfer, D., Leahy, M., Anctil, F. and Coulibaly, P., 2020. The Canadian Surface Prediction Archive (CaSPAr): A Platform to Enhance Environmental Modeling in Canada and Globally. Bulletin of the American Meteorological Society, 101(3), pp.E341-E356."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# This entire section is cookie-cutter template to allow calling the servers and instantiating the connection\n",
    "# to the WPS server. Do not modify this block.\n",
    "import datetime as dt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "from clisops.core import average, subset\n",
    "\n",
    "from ravenpy import Emulator\n",
    "from ravenpy.config import commands as rc\n",
    "from ravenpy.config.emulators import GR4JCN\n",
    "from ravenpy.extractors.forecasts import get_CASPAR_dataset\n",
    "from ravenpy.utilities import forecasting\n",
    "from ravenpy.utilities.testdata import get_file, open_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the warm-up file\n",
    "\n",
    "Here we tell the model that we want to forecast over the Salmon River catchment and provide its properties (area, lat/long, elevation). We will run it using the GR4JCN hydrological model and have provided some parameters. Other information on the forecast conditions is provided. Thr first step is to generate a hotstart file to prepare the model to generate forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the warmup period dates\n",
    "start_date_wu = dt.datetime(2010, 1, 1)\n",
    "end_date_wu = dt.datetime(2018, 6, 30)\n",
    "\n",
    "# Define the catchment contour. Here we use the Salmon River file we previously generated using the Delineator\n",
    "# in Tutorial Notebook 01.\n",
    "basin_contour = get_file(\"notebook_inputs/salmon_river.geojson\")\n",
    "\n",
    "# Define some of the catchment properties. Could also be replaced by a call to the properties WPS as in\n",
    "# the Tutorial Notebook 02.\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",
    "# Observed weather data for the Salmon river. We extracted this using Tutorial Notebook 03 and the\n",
    "# salmon_river.geojson file as the contour.\n",
    "ts = get_file(\"notebook_inputs/ERA5_weather_data_Salmon.nc\")\n",
    "\n",
    "# Set alternative names for netCDF variables\n",
    "alt_names = {\n",
    "    \"TEMP_MIN\": \"tmin\",\n",
    "    \"TEMP_MAX\": \"tmax\",\n",
    "    \"PRECIP\": \"pr\",\n",
    "}\n",
    "\n",
    "# Data types to extract from netCDF\n",
    "data_type = [\"TEMP_MAX\", \"TEMP_MIN\", \"PRECIP\"]\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"Latitude\": hru[\"latitude\"],\n",
    "        \"Longitude\": hru[\"longitude\"],\n",
    "    },\n",
    "}\n",
    "\n",
    "# Model configuration\n",
    "model_config_warmup = 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_wu,\n",
    "    EndDate=end_date_wu,\n",
    "    RunName=\"ESP_vs_NWP_warmup\",\n",
    ")\n",
    "\n",
    "# Run the model and get the outputs.\n",
    "out1 = Emulator(config=model_config_warmup).run()\n",
    "\n",
    "# Extract the path to the final states file that will be used as the next initial states\n",
    "hotstart = out1.files[\"solution\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dss = open_dataset(ts)\n",
    "dss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hindcasting using Climatological Ensemble Streamflow Prediction (ESP)\n",
    "Now that we have the hotstart file ready to go, we can configure our model for forecasting in climatology ESP mode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Date of the hindcast\n",
    "hdate = dt.datetime(2018, 7, 1)\n",
    "\n",
    "# Duration of the hindcast, in days\n",
    "duration = 7\n",
    "\n",
    "# Build a new model config:\n",
    "# Model configuration\n",
    "model_config_ESP = 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=hdate,\n",
    "    Duration=duration,\n",
    "    RunName=\"ESP_vs_NWP_ESPfcst\",\n",
    ")\n",
    "\n",
    "# Set the initial states of this new config to the correct values, i.e. the end of the previous forecast.\n",
    "model_config_ESP = model_config_ESP.set_solution(hotstart)\n",
    "\n",
    "# Simulate the climatological ESP:\n",
    "ESP_sims = forecasting.climatology_esp(config=model_config_ESP)\n",
    "\n",
    "# Show the results in an xarray dataset, ready to use:\n",
    "ESP_sims.hydrograph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now run the hindcast using Climatological ESP and retrieved the results. Let's take a look at the resulting forecast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Invent an observation so we can compute metrics later, and display as Qobs here. TODO: Add real streamflow data.\n",
    "qq = ESP_sims.hydrograph.q_sim[0, :, 0]\n",
    "\n",
    "# This is to be replaced with a call to the forecast graphing WPS as soon as merged.\n",
    "# model.q_sim.plot.line(\"b\", x=\"time\")\n",
    "ESP_sims.hydrograph.q_sim[:, :, 0].plot.line(\"b\", x=\"time\", add_legend=False)\n",
    "ESP_sims.hydrograph.q_sim[1, :, 0].plot.line(\"b\", x=\"time\", label=\"ESP forecasts\")\n",
    "qq.plot.line(\"r\", x=\"time\", label=\"observations\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hindcasting using archived weather forecasts from a weather forecast model\n",
    "\n",
    "In this next part, we will use the CaSPAr dataset (archived weather forecasts from Environment and Climate Change Canada) to forecast flows on the same period using the same hotstart file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the Forecast data from GEPS via CASPAR.\n",
    "# Take an extra day to ensure time-shift doesn't remove a part of our day\n",
    "ts_hindcast, _ = get_CASPAR_dataset(\"GEPS\", hdate - dt.timedelta(days=1))\n",
    "\n",
    "# Subset the data for the region of interest and take the mean to get a single vector\n",
    "with xr.set_options(keep_attrs=True):\n",
    "    ts_subset = subset.subset_shape(ts_hindcast, basin_contour).mean(\n",
    "        dim=(\"rlat\", \"rlon\")\n",
    "    )\n",
    "\n",
    "ts_subset = ts_subset.resample(time=\"6H\").nearest(\n",
    "    tolerance=\"1H\"\n",
    ")  # To make the timesteps identical accross the entire duration\n",
    "\n",
    "# We need to write the hindcast data as a file for Raven to be able to access it.\n",
    "fname = \"/tmp/hindcast.nc\"\n",
    "ts_subset.to_netcdf(fname)\n",
    "\n",
    "# We need to adjust the data_type and alt_names according to the data in the forecast:\n",
    "# Set alternative names for netCDF variables\n",
    "alt_names = {\n",
    "    \"TEMP_AVE\": \"tas\",\n",
    "    \"PRECIP\": \"pr\",\n",
    "}\n",
    "\n",
    "# Data types to extract from netCDF\n",
    "data_type = [\"TEMP_AVE\", \"PRECIP\"]\n",
    "\n",
    "# We will need to reuse this for GR4J. Update according to your needs. For example, here we will also pass\n",
    "# the catchment latitude and longitude as our CaSPAr data has been averaged at the catchment scale.\n",
    "# We also need to tell the model to deaccumulate the precipitation and shift it in time by 6 hours for our\n",
    "# catchment (UTC timezones):\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"Latitude\": hru[\"latitude\"],\n",
    "        \"Longitude\": hru[\"longitude\"],\n",
    "    },\n",
    "    \"PRECIP\": {\n",
    "        \"Deaccumulate\": True,\n",
    "        \"TimeShift\": -0.25,\n",
    "        \"LinearTransform\": {\n",
    "            \"scale\": 1000.0\n",
    "        },  # Since we are deaccumulating, we need to manually specify scale.\n",
    "    },  # Converting meters to mm (multiply by 1000).\n",
    "    \"TEMP_AVE\": {\n",
    "        \"TimeShift\": -0.25,\n",
    "    },\n",
    "}\n",
    "\n",
    "# Model configuration for forecasting, including correct start date and forecast duration\n",
    "model_config_fcst = GR4JCN(\n",
    "    params=[0.529, -3.396, 407.29, 1.072, 16.9, 0.947],\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            fname, data_type=data_type, alt_names=alt_names, data_kwds=data_kwds\n",
    "        )\n",
    "    ],\n",
    "    HRUs=[hru],\n",
    "    StartDate=hdate,\n",
    "    Duration=duration,\n",
    "    RunName=\"NB12_forecast_run\",\n",
    ")\n",
    "\n",
    "# Update the initial states\n",
    "model_config_fcst = model_config_fcst.set_solution(hotstart)\n",
    "\n",
    "# Generate the hindcast by providing all necessary information to generate virtual stations representing\n",
    "# the forecast members\n",
    "hindcast_sims = forecasting.hindcast_from_meteo_forecast(\n",
    "    model_config_fcst,\n",
    "    forecast=fname,\n",
    "    overwrite=True,\n",
    "    # We also need to provide the necessary information to create gauges inside the forecasting model:\n",
    "    data_kwds=data_kwds,\n",
    "    data_type=data_type,\n",
    "    alt_names=alt_names,\n",
    ")\n",
    "\n",
    "# Display the hydrographs\n",
    "display(hindcast_sims.hydrograph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hindcast_sims.hydrograph.q_sim[:, :, 0].plot.line(\"b\", x=\"time\", add_legend=False)\n",
    "hindcast_sims.hydrograph.q_sim[1, :, 0].plot.line(\"b\", x=\"time\", label=\"hindcasts\")\n",
    "qq.plot.line(\"r\", x=\"time\", label=\"observations\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model has run in forecast mode and we can now easily compare results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hindcast_sims.hydrograph.q_sim[:, :, 0].plot.line(\"b\", x=\"time\", add_legend=False)\n",
    "ESP_sims.hydrograph.q_sim[:, :, 0].plot.line(\"g\", x=\"time\", add_legend=False)\n",
    "ESP_sims.hydrograph.q_sim[1, :, 0].plot.line(\"g\", x=\"time\", label=\"ESP forecasts\")\n",
    "hindcast_sims.hydrograph.q_sim[1, :, 0].plot.line(\"b\", x=\"time\", label=\"hindcasts\")\n",
    "qq.plot.line(\"r\", x=\"time\", label=\"observations\")\n",
    "plt.legend(loc=\"upper left\")\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": 1
}