docs/notebooks/Running_HMETS_with_CANOPEX_dataset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forcing HMETS with the extended CANOPEX dataset\n",
    "\n",
    "Here we use ravenpy to launch the HMETS hydrological model and analyze the output. We also prepare and gather data directly from the CANOPEX dataset made available freely for all users."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cookie-cutter template necessary to provide the tools, packages and paths for the project. All notebooks\n",
    "# need this template (or a slightly adjusted one depending on the required packages)\n",
    "import datetime as dt\n",
    "import tempfile\n",
    "from pathlib import Path\n",
    "\n",
    "import pandas as pd\n",
    "import spotpy\n",
    "import xarray as xr\n",
    "\n",
    "from ravenpy.config import commands as rc\n",
    "from ravenpy.config.emulators import HMETS\n",
    "from ravenpy.utilities.calibration import SpotSetup\n",
    "from ravenpy.utilities.testdata import get_file\n",
    "\n",
    "# Make a temporary folder\n",
    "tmp = Path(tempfile.mkdtemp())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DATA MAIN SOURCE - DAP link to CANOPEX dataset.\n",
    "CANOPEX_DAP = \"https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/birdhouse/ets/Watersheds_5797_cfcompliant.nc\"\n",
    "ds = xr.open_dataset(CANOPEX_DAP)\n",
    "\n",
    "# Explore the dataset:\n",
    "display(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We could explore the dataset and find a watershed of interest, but for now, let's pick one at random\n",
    "# from the dataset:\n",
    "watershedID = 5600\n",
    "\n",
    "# And show what it includes:\n",
    "ts = ds.isel({\"watershed\": watershedID})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's write the file to disk to make it more efficient to retrieve:\n",
    "fname = tmp / \"CANOPEX_extracted.nc\"\n",
    "ts.to_netcdf(fname)\n",
    "ds.close()\n",
    "ts.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# With this info, we can gather some properties from the CANOPEX database. This same database is used for\n",
    "# regionalization, so let's query it there where more information is available:\n",
    "tmp = pd.read_csv(get_file(\"regionalisation_data/gauged_catchment_properties.csv\"))\n",
    "\n",
    "basin_area = float(tmp[\"area\"][watershedID])\n",
    "basin_latitude = float(tmp[\"latitude\"][watershedID])\n",
    "basin_longitude = float(tmp[\"longitude\"][watershedID])\n",
    "basin_elevation = float(tmp[\"elevation\"][watershedID])\n",
    "basin_name = ds.watershed.data\n",
    "\n",
    "print(\"Basin name: \", basin_name)\n",
    "print(\"Latitude: \", basin_latitude, \" °N\")\n",
    "print(\"Area: \", basin_area, \" km^2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we might have the model and data, but we don't have model parameters! We need to calibrate. This next snippets show how to configure the model and the calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will also calibrate on only a subset of the years for now to keep the computations faster in this notebook.\n",
    "start_calib = dt.datetime(1998, 1, 1)\n",
    "end_calib = dt.datetime(1999, 12, 31)\n",
    "\n",
    "# General parameters depending on the data source. We can find them by exploring the CANOPEX dataset in the\n",
    "# cells above.\n",
    "data_type = [\"TEMP_MAX\", \"TEMP_MIN\", \"PRECIP\"]\n",
    "\n",
    "alt_names = {\n",
    "    \"TEMP_MIN\": \"tasmin\",\n",
    "    \"TEMP_MAX\": \"tasmax\",\n",
    "    \"PRECIP\": \"pr\",\n",
    "}\n",
    "\n",
    "hru = {}\n",
    "hru = dict(\n",
    "    area=basin_area,\n",
    "    elevation=basin_elevation,\n",
    "    latitude=basin_latitude,\n",
    "    longitude=basin_longitude,\n",
    "    hru_type=\"land\",\n",
    ")\n",
    "\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"latitude\": hru[\"latitude\"],\n",
    "        \"longitude\": hru[\"longitude\"],\n",
    "    }\n",
    "}\n",
    "# Set the evaluation metrics to be calculated by Raven\n",
    "eval_metrics = (\"NASH_SUTCLIFFE\",)\n",
    "\n",
    "model_config = HMETS(\n",
    "    ObservationData=[\n",
    "        rc.ObservationData.from_nc(fname, alt_names=\"discharge\", station_idx=1)\n",
    "    ],\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            fname,\n",
    "            station_idx=1,\n",
    "            data_type=data_type,  # Note that this is the list of all the variables\n",
    "            alt_names=alt_names,  # Note that all variables here are mapped to their names in the netcdf file.\n",
    "            data_kwds=data_kwds,\n",
    "        )\n",
    "    ],\n",
    "    HRUs=[hru],\n",
    "    StartDate=start_calib,\n",
    "    EndDate=end_calib,\n",
    "    RunName=\"CANOPEX_test\",\n",
    "    EvaluationMetrics=eval_metrics,\n",
    "    RainSnowFraction=\"RAINSNOW_DINGMAN\",\n",
    "    SuppressOutput=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The model parameters bounds can either be set independently or we can use the defaults.\n",
    "low_params = (\n",
    "    0.3,\n",
    "    0.01,\n",
    "    0.5,\n",
    "    0.15,\n",
    "    0.0,\n",
    "    0.0,\n",
    "    -2.0,\n",
    "    0.01,\n",
    "    0.0,\n",
    "    0.01,\n",
    "    0.005,\n",
    "    -5.0,\n",
    "    0.0,\n",
    "    0.0,\n",
    "    0.0,\n",
    "    0.0,\n",
    "    0.00001,\n",
    "    0.0,\n",
    "    0.00001,\n",
    "    0.0,\n",
    "    0.0,\n",
    ")\n",
    "high_params = (\n",
    "    20.0,\n",
    "    5.0,\n",
    "    13.0,\n",
    "    1.5,\n",
    "    20.0,\n",
    "    20.0,\n",
    "    3.0,\n",
    "    0.2,\n",
    "    0.1,\n",
    "    0.3,\n",
    "    0.1,\n",
    "    2.0,\n",
    "    5.0,\n",
    "    1.0,\n",
    "    3.0,\n",
    "    1.0,\n",
    "    0.02,\n",
    "    0.1,\n",
    "    0.01,\n",
    "    0.5,\n",
    "    2.0,\n",
    ")\n",
    "\n",
    "# Setup the spotpy optimizer\n",
    "spot_setup = SpotSetup(\n",
    "    config=model_config,\n",
    "    low=low_params,\n",
    "    high=high_params,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can run the optimizer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll definitely want to adjust the random seed and number of model evaluations:\n",
    "model_evaluations = (\n",
    "    50  # This is to keep computing time fast for the demo, increase as necessary\n",
    ")\n",
    "\n",
    "# Setup the spotpy sampler with the method, the setup configuration, a run name and other options. Please refer to\n",
    "# the spotpy documentation for more options. We recommend sticking to this format for efficiency of most applications.\n",
    "sampler = spotpy.algorithms.dds(\n",
    "    spot_setup,\n",
    "    dbname=\"CANOPEX_test\",\n",
    "    dbformat=\"ram\",\n",
    "    save_sim=False,\n",
    ")\n",
    "\n",
    "# Launch the actual optimization. Multiple trials can be launched, where the entire process is repeated and\n",
    "# the best overall value from all trials is returned.\n",
    "sampler.sample(model_evaluations, trials=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the model diagnostics\n",
    "diag = spot_setup.diagnostics\n",
    "\n",
    "# Print the NSE and the parameter set in 2 different ways:\n",
    "print(\"Nash-Sutcliffe value is: \" + str(diag[\"DIAG_NASH_SUTCLIFFE\"]))\n",
    "\n",
    "# Get all the values of each iteration\n",
    "results = sampler.getdata()\n",
    "\n",
    "# Get the raw resutlts directly in an array\n",
    "params = spotpy.analyser.get_best_parameterset(results)[0]\n",
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this stage, we have calibrated the model on the observations for the desired dates. Now, let's run the model on a longer time period and look at the hydrograph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ravenpy import Emulator\n",
    "\n",
    "conf = model_config.set_params(params)\n",
    "conf.suppress_output = False\n",
    "out = Emulator(conf).run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `hydrograph` and `storage` outputs are netCDF files storing the time series. These files are opened by default using `xarray`, which provides convenient and powerful time series analysis and plotting tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = out.hydrograph.q_sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can also get statistics from the data directly.\n",
    "print(\"Max: \", q.max().values)\n",
    "print(\"Mean: \", q.mean().values)\n",
    "print(\n",
    "    \"Monthly means: \",\n",
    "    q.groupby(\"time.month\").mean(dim=\"time\").values,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot the simulated hydrograph\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "\n",
    "register_matplotlib_converters()\n",
    "q.plot()"
   ]
  }
 ],
 "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": 4
}