{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(_selle_morvillers)=\n",
    "# Modelling groundwater levels and river flows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is taken from the Gardenia tutorial {cite:p}`2013:thiery_didacticiel`.\n",
    "\n",
    "This example simulates river flows and groundwater levels in a single watershed with \n",
    "daily observed river flows and monthly observed groundwater levels. The following data\n",
    "are available:\n",
    "\n",
    "- Observed daily river flow of the Selle river at Plachy from October 1989 to July 2003 (m3/s)\n",
    "- Observed monthly groundwater levels at Morvillers from September 1985 to July 2003 (mNGF)\n",
    "- Daily rainfall and potential evapotranspiration from August 1985 to July 2003 (mm/day)\n",
    "\n",
    "The Selle River basin is a sub-basin of the Somme River basin. Its drainage area is equal to 524 km².\n",
    "\n",
    "A simulation starts from 1989-01-01 until 2003-07-31, date of the rainfall data\n",
    "last record. A warmup period from 1985-01-01 to 1988-12-32 repeated three times\n",
    "initialize the model:\n",
    "\n",
    "Four parameters will be optimized through bound-constrained optimization: soil\n",
    "water capacity, runoff/seepage partition coefficient, transfer halflife time,\n",
    "and baseflow halflife time.\n",
    "\n",
    "Two will be optimized through regression: base level and groundwater storage coefficient.\n",
    "\n",
    "\n",
    "Optimization of parameters uses 200 iterations at maximum and the Nash-Sutcliff\n",
    "coefficient as river flow metric in the objective function. Metrics are computed\n",
    "from 1989-01-01 to the end of the simulation (2003-07-31).\n",
    "\n",
    "In the objective function, the weight for the river flow metric is equal to 5\n",
    "and the weight for the groundwater level metric (Nash-Sutcliff coefficient) is\n",
    "equal to 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TOML configuration file\n",
    "\n",
    "```{literalinclude} model.toml\n",
    ":language: toml\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import rameau as rm\n",
    "\n",
    "# Load model from a toml file\n",
    "model = rm.Model.from_toml(f\"model.toml\")\n",
    "\n",
    "# Run a simulation with parameters defined in the toml file\n",
    "sim = model.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the observed and simulated river flow time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Get the riverflow simulation\n",
    "riv_sim = sim.get_output(\"riverflow\")\n",
    "\n",
    "# Get the riverflow observation\n",
    "riv_obs = model.get_input(\"riverobs\")\n",
    "\n",
    "# Plot river flow\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":riv_obs.iloc[:, 0],\n",
    "        \"Sim\":riv_sim.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1989-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed riverflow (Selle à Plachy)\",\n",
    "    ylabel=\"m3/s\",\n",
    "    style=[\"-\", \"+\"]\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the observed and groundwater level time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the watertable simulation\n",
    "wtl_sim = sim.get_output(\"watertable\")\n",
    "\n",
    "# Get the watertable observation\n",
    "wtl_obs = model.get_input(\"groundwaterobs\")\n",
    "\n",
    "# Plot watertable\n",
    "df = pd.DataFrame({\"Sim\":wtl_sim.iloc[:, 0], \"Obs\":wtl_obs.iloc[:, 0]})\n",
    "df.loc[\"1989-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed watertable level (Selle à Plachy)\",\n",
    "    ylabel=\"m NGF\",\n",
    "    style=['-', '+']\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the riverflow metrics\n",
    "scores = sim.get_metrics(\"riverflow\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the watertable metrics\n",
    "scores = sim.get_metrics(\"watertable\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the optimization settings already that have been loaded through the toml file,\n",
    "# in the [optimization] section.\n",
    "for key, value in model.optimization_settings.items():\n",
    "    print(f'{key} = {value}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a first optimization with the already loaded settings\n",
    "sim_opt1 = model.run_optimization()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the riverflow simulation from optimization 1 \n",
    "riv_sim_opt1 = sim_opt1.get_output(\"riverflow\")\n",
    "\n",
    "# Plot river flow\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":riv_obs.iloc[:, 0],\n",
    "        \"Sim\":riv_sim.iloc[:, 0],\n",
    "        \"Opt\":riv_sim_opt1.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1989-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed riverflow (Selle à Plachy)\",\n",
    "    ylabel=\"m3/s\",\n",
    "    style=[\"+\", \"-\", \"-\"]\n",
    ")\n",
    "plt.show() \n",
    "\n",
    "# Print the riverflow optimization metrics\n",
    "scores = sim_opt1.get_opti_metrics(\"riverflow\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the watertable simulation from optimization 1 \n",
    "wtl_sim_opt1 = sim_opt1.get_output(\"watertable\")\n",
    "\n",
    "# Plot watertable\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":wtl_obs.iloc[:, 0],\n",
    "        \"Sim\":wtl_sim.iloc[:, 0],\n",
    "        \"Opt\":wtl_sim_opt1.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.loc[\"1989-01-01\":, :].plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed watertable level (Selle à Plachy)\",\n",
    "    ylabel=\"m NGF\",\n",
    "    style=[\"+\", \"-\", \"-\"]\n",
    ")\n",
    "plt.show() \n",
    "\n",
    "# Print the groundwater level optimization metrics\n",
    "scores = sim_opt1.get_opti_metrics(\"watertable\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now run a second optimization but with different options\n",
    "# (Gives almost the same results at the end...)\n",
    "sim_opt2 = model.run_optimization(\n",
    "    maxit=400,\n",
    "    river_objective_function=\"kge_2012\",\n",
    "    transformation=\"square root\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the riverflow simulation from optimization 2 \n",
    "riv_sim_opt2 = sim_opt1.get_output(\"riverflow\")\n",
    "\n",
    "# Plot river flow\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"Obs\":riv_obs.iloc[:, 0],\n",
    "        \"Sim\":riv_sim.iloc[:, 0],\n",
    "        \"Opt\":riv_sim_opt1.iloc[:, 0],\n",
    "        \"Opt 2\":riv_sim_opt2.iloc[:, 0],\n",
    "    }\n",
    ")\n",
    "df.plot(\n",
    "    grid=True,\n",
    "    title=\"Simulated and observed riverflow (Selle à Plachy)\",\n",
    "    ylabel=\"m3/s\",\n",
    "    style=[\"+\", \"-\", \"-\"]\n",
    ")\n",
    "plt.show() \n",
    "\n",
    "# Print the riverflow optimization metrics for optimization 2\n",
    "scores = sim_opt2.get_opti_metrics(\"riverflow\")\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update the parameters of the model\n",
    "model.tree = sim_opt2.tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print default options for forecast\n",
    "for key, value in model.forecast_settings.items():\n",
    "    print(f'{key} = {value}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a forecast with options\n",
    "import datetime\n",
    "\n",
    "sim_for = model.run_forecast(\n",
    "    scope=datetime.timedelta(days=90),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now print options for forecast that have been defined for the simulation sim_for\n",
    "# The default emission date corresponds to the last date of the known meteorological data\n",
    "for key, value in sim_for.forecast_settings.items():\n",
    "    print(f'{key} = {value}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get historical river flow\n",
    "hist = sim_for.get_output(\"riverflow\")\n",
    "\n",
    "# Get forecast river flow\n",
    "fore = sim_for.get_forecast_output(\"riverflow\")\n",
    "\n",
    "# Plot the spaghettis\n",
    "riv_obs.columns = [\"Obs\"]\n",
    "s = datetime.datetime(2003, 1, 1)\n",
    "e = hist.index[-1] + sim_for.forecast_settings.scope\n",
    "fig, ax = plt.subplots()\n",
    "riv_obs.loc[s:e,:].plot(ax=ax)\n",
    "hist.loc[s:e, :].plot(ax=ax, legend=False)\n",
    "fore.loc[s:e, :].plot(ax=ax, legend=False)\n",
    "ax.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a forecast with other options\n",
    "sim_for = model.run_forecast(\n",
    "    scope=datetime.timedelta(days=90),\n",
    "    quantiles_output=True,\n",
    "    norain=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get historical river flow\n",
    "hist = sim_for.get_output(\"riverflow\")\n",
    "\n",
    "# Get forecast river flow\n",
    "fore = sim_for.get_forecast_output(\"riverflow\")\n",
    "fore = fore.xs(key=\"Selle at Plachy\", level=\"watersheds\", axis=1)\n",
    "\n",
    "# Plot the spaghettis\n",
    "riv_obs.columns = [\"Obs\"]\n",
    "hist.columns = [\"Hist\"]\n",
    "s = datetime.datetime(2003, 1, 1)\n",
    "e = hist.index[-1] + sim_for.forecast_settings.scope\n",
    "fig, ax = plt.subplots()\n",
    "riv_obs.loc[s:e,:].plot(ax=ax)\n",
    "hist.loc[s:e, :].plot(ax=ax)\n",
    "fore.loc[s:e, :].plot(ax=ax)\n",
    "ax.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run a forecast with halflife correction.\n",
    "model.tree.watersheds[0].forecast_correction.river.halflife = 30\n",
    "sim_hl = model.run_forecast(\n",
    "    scope=datetime.timedelta(days=90),\n",
    "    quantiles_output=True,\n",
    "    norain=True,\n",
    "    correction=\"halflife\"\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get historical river flow\n",
    "hist = sim_hl.get_output(\"riverflow\")\n",
    "\n",
    "# Get forecast river flow\n",
    "fore = sim_hl.get_forecast_output(\"riverflow\")\n",
    "fore = fore.xs(key=\"Selle at Plachy\", level=\"watersheds\", axis=1)\n",
    "\n",
    "# Plot the spaghettis\n",
    "riv_obs.columns = [\"Obs\"]\n",
    "hist.columns = [\"Hist\"]\n",
    "s = datetime.datetime(2003, 1, 1)\n",
    "e = hist.index[-1] + sim_for.forecast_settings.scope\n",
    "fig, ax = plt.subplots()\n",
    "riv_obs.loc[s:e,:].plot(ax=ax)\n",
    "hist.loc[s:e, :].plot(ax=ax)\n",
    "fore.loc[s:e, :].plot(ax=ax)\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev",
   "language": "python",
   "name": "python3"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}