docs/notebooks/05_Advanced_RavenPy_configuration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 05 - Advanced RavenPy configuration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will explore alternative ways to setup a Raven model and how to parameterize and customize a raven-based hydrological model\n",
    "\n",
    "## Running Raven using pre-existing configuration files\n",
    "\n",
    "To run Raven, we need configuration (`.rvX`) files defining hydrological processes, watersheds and meteorological data. If you already have those configuration files ready, or want to see how to import an existing Raven model into PAVICS-Hydro, this tutorial is for you. It shows how to run Raven from a Python programming environment using [RavenPy](https://ravenpy.readthedocs.io/en/latest/).\n",
    "\n",
    "Let's start by importing some utilities that will make our life easier to get data on the servers. If you already have raven model setups, you could simply upload the files here and create your own \"config\" list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility that simplifies getting data hosted on the remote PAVICS-Hydro data server.\n",
    "from ravenpy.utilities.testdata import get_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A note on datasets\n",
    "\n",
    "For this part of the tutorial, we will use pre-existing datasets that are hosted on the PAVICS-Hydro servers to setup the Raven model. This means that the .rv files are all built and the forcing file already exists. We could apply all of the same logic to a RavenPy model we would have built at the previous step, but this way lets us show that we can also work on an imported model. Let's import the configuration files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the .rv files. It could also be the .rv files returned from the previous notebook, but here we are using a new basin that contains observed streamflow\n",
    "# to make the calibration possible in the next notebook. Note that these configuration files also include links to the\n",
    "# required hydrometeorological database (NetCDF file).\n",
    "config = [\n",
    "    get_file(f\"raven-gr4j-cemaneige/raven-gr4j-salmon.{ext}\")\n",
    "    for ext in [\"rvt\", \"rvc\", \"rvi\", \"rvh\", \"rvp\"]\n",
    "]\n",
    "config"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So \"config\" is just a set of paths to the various .rvX files (.rvt, .rvc, .rvi. .rvh and .rvp). Therefore, if you have your own .rv files that describe your model, you can upload them and replace \"config\" with your own files!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a hydrological model on-the-fly using existing configuration files.\n",
    "\n",
    "Here we create a Raven model instance, configuring it using the pre-defined configuration files and running it by providing the full path to the NetCDF driving datasets. The configuration we provide is for a GR4J-CN model emulator that Raven will run for us. We provide the configuration files for GR4J-CN as well as the forcing data (precipitation, temperature, observed streamflow, etc.) that will be used to run the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ravenpy import OutputReader\n",
    "from ravenpy.ravenpy import run\n",
    "\n",
    "run_name = \"raven-gr4j-salmon\"  # As can be seen in the config above, this is the name of the .rvX files.\n",
    "configdir = config[\n",
    "    0\n",
    "].parent  # We can get the path to the folder containing the .rvX files this way\n",
    "\n",
    "# Run the model and get the path to outputs\n",
    "outputs_path = run(modelname=run_name, configdir=configdir, overwrite=True)\n",
    "\n",
    "# Note. The modelname parameter can be confusing. You need to give the FILES extension name (run_name in our case),\n",
    "# not the name of the model.\n",
    "\n",
    "outputs_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the output files at the output_path\n",
    "\n",
    "outputs = OutputReader(run_name=None, path=outputs_path)  # Get the outputs\n",
    "# Note. We set up the run_name to None, because we didn't rename the output files. If you gave a different name to your file\n",
    "# compared to the one above, you should change the run_name value to this new name. It's important though that you keep the end\n",
    "# of the filename the same\n",
    "\n",
    "# Show the list of files that were retrived by the OutputReader\n",
    "outputs.files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model should have run! But you also might have seen some warnings that Raven is giving us, depending on the input files used:\n",
    "\n",
    "- Some might be saying that we are providing rain and snow independently, but in the configuration files, we are asking the model to recompute the separation using an algorithm based on total precipitation and air temperature. This is OK, and we can live with this (alternatively, we could reconfigure the model to remove this but that will be for another notebook!).\n",
    "\n",
    "- Others could be saying that we supply PET data, but the model is configured to compute PET from the available temperature and latitude/longitude data. This is also acceptable to us for now, so these warnings can be disregarded.\n",
    "\n",
    "- And others might simply explain that our configuration provided some parameters but others were computed internally based on our parameter set rather than being explicitly set in our configuration, which is OK.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluating the model response\n",
    "\n",
    "That's it! The code above has launched the GR4J-CN model using weather data and the configuration we provided. There are many other options we could provide, but for now we left everything to the default options to keep things simple. We will explore those in a future tutorial as well.\n",
    "\n",
    "Now, let's look at the modeled hydrographs. Note that there is a \"q_obs\" hydrograph, representing the observations we provided ourselves. This is to facilitate the comparison between observations and simulations, and it is not required per se to run the model. The \"q_sim\" variable is the simulated streamflow and is the one we are interested in.\n",
    "\n",
    "Note that RavenPy assumes that model outputs are always saved in netCDF format, and relies on [xarray](http://xarray.pydata.org/en/stable/) to access data.\n",
    "\n",
    "To see results, we must first tell the model to read them from the files Raven has written in the output folder:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the simulated streamflow using xarray's built-in plotting tool, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outputs.hydrograph.q_sim.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also now have access to diagnostics! This is because along with the simulated discharge, the model has access to observed discharge to compute error metrics such as RMSE and NSE. Let's see where the file has been generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"-----------------DIAGNOSTICS-----------------\")\n",
    "print(outputs.diagnostics)\n",
    "print(\"\")\n",
    "\n",
    "print(\"-----------------NASH_SUTCLIFFE-----------------\")\n",
    "print(outputs.diagnostics[\"DIAG_NASH_SUTCLIFFE\"])\n",
    "print(\"\")\n",
    "\n",
    "print(\"-----------------RMSE-----------------\")\n",
    "print(outputs.diagnostics[\"DIAG_RMSE\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the Nash-Sutcliffe value is quite poor. This is due to the short simulation period in the configuration (see the hydrograph above!) and the lack of a spin-up period, combined to a poor parameter set choice. We will improve upon all of these shortcomings in the next notebooks!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced RavenPy configuration options\n",
    "\n",
    "Raven can perform many operations and has multiple configuration options. Here we provide a list of configuration options to explore which you can eventually use to taylor the codes to your own specifications. These can only be run on RavenPy-built hydrological models, and will not operate on Raven models imported by users since those configuration files are not modifiable for the time being.\n",
    "\n",
    "We will give an overview of the various configuration keywords after this code block, but users should read the Raven documentation for more options for each of these processes.\n",
    "\n",
    "Let's first define some variables we will need for all of our tests:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get required packages\n",
    "import datetime as dt\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from ravenpy import Emulator\n",
    "from ravenpy.config import commands as rc\n",
    "from ravenpy.config.emulators import GR4JCN\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 alternate variable names in the timeseries data file\n",
    "alt_names = {\n",
    "    \"TEMP_MIN\": \"tmin\",\n",
    "    \"TEMP_MAX\": \"tmax\",\n",
    "    \"PRECIP\": \"pr\",\n",
    "}\n",
    "\n",
    "# Provide the type of data made available to Raven\n",
    "data_type = [\"TEMP_MAX\", \"TEMP_MIN\", \"PRECIP\"]\n",
    "\n",
    "# Prepare the catchment properties\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",
    "# Add some information regarding station data\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"latitude\": hru[\"latitude\"],\n",
    "        \"longitude\": hru[\"longitude\"],\n",
    "    }\n",
    "}\n",
    "\n",
    "# Start and end dates of the simulation\n",
    "start_date = dt.datetime(1985, 1, 1)\n",
    "end_date = dt.datetime(1990, 1, 1)\n",
    "\n",
    "# Set parameters\n",
    "parameters = [0.529, -3.396, 407.29, 1.072, 16.9, 0.947]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now perform a \"basic\" run, with no modifications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the model (See Notebook 04 for more details on implementation)\n",
    "m = GR4JCN(\n",
    "    params=parameters,\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            ts,\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_date,\n",
    "    EndDate=end_date,\n",
    "    RunName=\"NB05_test1\",\n",
    "    # GlobalParameter={\"AVG_ANNUAL_RUNOFF\": 208.480},\n",
    ")\n",
    "\n",
    "# Run the model and get the outputs.\n",
    "outputs1 = Emulator(m).run()\n",
    "\n",
    "# Plot the generated hydrograph\n",
    "outputs1.hydrograph.q_sim.plot.line(x=\"time\", label=\"Base case\")\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now run another model by adding some other properties. To start, we can add some Global Parameters to the model to make Raven adjust the simulations based on the information we provide. Some options of Global Parameters are indicated here, but more can be found in the official Raven documentation.\n",
    "\n",
    "Examples of GlobalParameter options (Note that some are only available for certain models and others can be mutually exclusive. Please refer to the documentation for this type of adjustment):\n",
    "\n",
    "### Temperature interval of transformation between rain and snow. Set the midpoint of the range and the width of the range, in degrees C:\n",
    "\"RAINSNOW_TEMP\": midpoint_temp  // Ex: \"RAINSNOW_TEMP\": -1.0\n",
    "\n",
    "\"RAINSNOW_DELTA\": delta_temp    // Ex: \"RAINSNOW_DELTA\": 3.0\n",
    "\n",
    "### Maximum liquid water content of snow, as a percentage of SWE (0-1). Usually ~0.05.\n",
    "\"SNOW_SWI\": saturation // Ex: \"SNOW_SWI\": 0.1\n",
    "\n",
    "### Average annual snow for the entire watershed in mm of SWE. Used in CemaNeige.\n",
    "\"AVG_ANNUAL_SNOW\": average_snow_per_year // Ex: \"AVG_ANNUAL_SNOW\": 400.0\n",
    "\n",
    "There are many others, but this should clarify the implementation. Let's try some of them out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the model (See Notebook 04 for more details on implementation)\n",
    "m = GR4JCN(\n",
    "    params=parameters,\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            ts,\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_date,\n",
    "    EndDate=end_date,\n",
    "    RunName=\"NB05_test2\",\n",
    "    GlobalParameter={\"AVG_ANNUAL_SNOW\": 350.0},\n",
    ")\n",
    "\n",
    "# Run the model and get the outputs.\n",
    "outputs2 = Emulator(m).run()\n",
    "\n",
    "# Plot the generated hydrograph\n",
    "outputs1.hydrograph.q_sim.plot.line(x=\"time\", label=\"Base case\")\n",
    "outputs2.hydrograph.q_sim.plot.line(x=\"time\", label=\"With AVG_ANNUAL_SNOW\")\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also adjust the time series data to play with the scaling of units.\n",
    "\n",
    "By default, RavenPy and Raven will detect units from the forcing data netcdf files. However, in some instances, units might be lacking, or their format might require some tinkering. One such case is for precipitation data that is cumulative in the netcdf file. In these cases, Raven can decumulate the precipitation, but the scaling might lead to undesirable results. For this reason, it is highly recommended to pass the scaling and offsetting variables directly. To do so, add some context in the data_kwds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add some information regarding station data\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"Latitude\": hru[\"latitude\"],\n",
    "        \"Longitude\": hru[\"longitude\"],\n",
    "    },\n",
    "    # HOW TO PROCESS THE PRECIPITATION DATA: For the Precip variable, we tell Raven we want to Deaccumulate\n",
    "    # values, shift them in time by 6 hours (for UTC time zone management), and then apply a linear transform\n",
    "    # to the values to get new scaled values. The linear transform can take two inputs:\n",
    "    #     \"scale\"   is the \"a\" variable in the linear relationship y = ax + b. Usually used to multiply precipitation.\n",
    "    #     \"offset\"  is the \"b\" variable in the linear relationship y = ax + b. Usually used to convert temperatures(K to °C)\n",
    "    \"PRECIP\": {\n",
    "        \"Deaccumulate\": True,\n",
    "        \"TimeShift\": -0.25,\n",
    "        \"LinearTransform\": {\n",
    "            \"scale\": 1000.0  # # Converting meters to mm (multiply by 1000).\n",
    "        },\n",
    "    },\n",
    "    \"TEMP_AVE\": {\n",
    "        \"TimeShift\": -0.25,\n",
    "    },\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our example, our precipitation is not actually accumulated and the timestep is daily, so we don't need the \"Deaccumulate\" or the \"TimeShift\" parameters. So let's generate a new data_kwds that is applicable in our case. More complex cases that require \"Deaccumulate\" and \"TimeShift\" will be presented in later notebooks that use accumulated precipitation in forecasting applications, in Notebook 12."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add some information regarding station data\n",
    "data_kwds = {\n",
    "    \"ALL\": {\n",
    "        \"elevation\": hru[\"elevation\"],\n",
    "        \"Latitude\": hru[\"latitude\"],\n",
    "        \"Longitude\": hru[\"longitude\"],\n",
    "    },\n",
    "    # Let's simulate a very rough estimation of the impacts of climate change where precipitation is expected\n",
    "    # to increase by 10% and temperatures to increase by 3°C. This will be applied to all data on the entire\n",
    "    # period and is thus not realistic. We will explore more realistic methods in Notebook 08.\n",
    "    \"PRECIP\": {\"LinearTransform\": {\"scale\": 1.1}},\n",
    "    \"TEMP_AVE\": {\"LinearTransform\": {\"offset\": 3.0}},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use this new setup to generate another series of streamflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the model (See Notebook 04 for more details on implementation)\n",
    "m = GR4JCN(\n",
    "    params=parameters,\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            ts,\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_date,\n",
    "    EndDate=end_date,\n",
    "    RunName=\"NB05_test3\",\n",
    "    GlobalParameter={\"AVG_ANNUAL_SNOW\": 350.0},\n",
    ")\n",
    "\n",
    "# Run the model and get the outputs.\n",
    "outputs3 = Emulator(m).run()\n",
    "\n",
    "# Plot the generated hydrograph\n",
    "outputs1.hydrograph.q_sim.plot.line(x=\"time\", label=\"Base case\")\n",
    "outputs2.hydrograph.q_sim.plot.line(x=\"time\", label=\"With AVG_ANNUAL_SNOW\")\n",
    "outputs3.hydrograph.q_sim.plot.line(x=\"time\", label=\"With AVG_ANNUAL_SNOW and Scaling\")\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the scaling increased the flows almost everywhere except in the first year which is the warm-up period.\n",
    "\n",
    "Other options that can be implemented are indicated here, although more exist and are documented in the official Raven manual.\n",
    "\n",
    "\n",
    "\n",
    "### RainSnowFraction:\n",
    "Algorithm to use to separate the total precipitation into rainfall and snowfall.\n",
    "\n",
    "Ex: RainSnowFraction='RAINSNOW_DINGMAN'\n",
    "\n",
    "### Evaporation\n",
    "Evaporation: Formula to use to compute the evapotranspiration from the land HRUs.\n",
    "\n",
    "Ex: Evaporation=\"PET_OUDIN\"\n",
    "\n",
    "### Suppress model outputs / files\n",
    "Boolean that indicates if you wish for Raven to provide information after the model evaluation by writing to file. For a single run this can be left to **False**, but for calibration and other intensive tasks, it is faster to leave it to **True**.\n",
    "\n",
    "Ex: SuppressOutputs=True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's see how to implement these commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the model (See Notebook 04 for more details on implementation)\n",
    "m = GR4JCN(\n",
    "    params=parameters,\n",
    "    Gauge=[\n",
    "        rc.Gauge.from_nc(\n",
    "            ts,\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_date,\n",
    "    EndDate=end_date,\n",
    "    RunName=\"NB05_test3\",\n",
    "    GlobalParameter={\"AVG_ANNUAL_SNOW\": 350.0},\n",
    "    RainSnowFraction=\"RAINSNOW_DINGMAN\",\n",
    "    Evaporation=\"PET_HARGREAVES_1985\",\n",
    "    SuppressOutput=False,  # We can't read the hydrographs if they are not written to disk, so set to False here.\n",
    ")\n",
    "\n",
    "# Run the model and get the outputs.\n",
    "outputs4 = Emulator(m).run()\n",
    "\n",
    "# Plot the generated hydrograph\n",
    "outputs1.hydrograph.q_sim.plot.line(x=\"time\", label=\"Base case\")\n",
    "outputs2.hydrograph.q_sim.plot.line(x=\"time\", label=\"With AVG_ANNUAL_SNOW\")\n",
    "outputs3.hydrograph.q_sim.plot.line(x=\"time\", label=\"With AVG_ANNUAL_SNOW and Scaling\")\n",
    "outputs4.hydrograph.q_sim.plot.line(\n",
    "    x=\"time\", label=\"With AVG_ANNUAL_SNOW, Scaling and Options\"\n",
    ")\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A note on the above results\n",
    "\n",
    "We can see that the results change significantly according to the options we have passed, namely the evaporation algorithm modified the hydrograph quite significantly. However, this is caused by the fact that the parameter set we have used has not been calibrated using this PET method, and thereore the model cannot be expected to perform as well. This means that when using these model options, it is important to recalibrate the model parameters such that they represent the actual model being used!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Finally, we can also ask Raven to supply custom outputs using this line in the model configuration:\n",
    "\n",
    "CustomOutput=rc.CustomOutput() and by providing a list of desired pre-processed variables. Here we ask for the yearly average of precipitation over the entire watershed:\n",
    "\n",
    "CustomOutput=rc.CustomOutput(\"YEARLY\", \"AVERAGE\", \"PRECIP\", \"ENTIRE_WATERSHED\")\n",
    "\n",
    "Please see the documentation for more details on using custom outputs.\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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}