{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Basic Model\n",
    "\n",
    "In this example application it is shown how a simple time series model can be developed to simulate groundwater levels. The recharge (calculated as precipitation minus evaporation) is used as the explanatory time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "import pastas as ps\n",
    "\n",
    "ps.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Importing the dependent time series data\n",
    "In this codeblock a time series of groundwater levels is imported using the `read_csv` function of `pandas`. As `pastas` expects a `pandas` `Series` object, the data is squeezed. To check if you have the correct data type (a `pandas Series` object), you can use `type(oseries)` as shown below. \n",
    "\n",
    "The following characteristics are important when importing and preparing the observed time series:\n",
    "- The observed time series are stored as a `pandas Series` object.\n",
    "- The time step can be irregular."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import groundwater time seriesm and squeeze to Series object\n",
    "gwdata = pd.read_csv(\n",
    "    \"data/head_nb1.csv\", parse_dates=[\"date\"], index_col=\"date\"\n",
    ").squeeze()\n",
    "print(\"The data type of the oseries is: %s\" % type(gwdata))\n",
    "\n",
    "# Plot the observed groundwater levels\n",
    "gwdata.plot(style=\".\", figsize=(10, 4))\n",
    "plt.ylabel(\"Head [m]\")\n",
    "plt.xlabel(\"Time [years]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import the independent time series\n",
    "Two explanatory series are used: the precipitation and the potential evaporation. These need to be `pandas Series` objects, as for the observed heads.\n",
    "\n",
    "Important characteristics of these time series are:\n",
    "- All series are stored as `pandas Series` objects.\n",
    "- The series may have irregular time intervals, but then it will be converted to regular time intervals when creating the time series model later on.\n",
    "- It is preferred to use the same length units as for the observed heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import observed precipitation series\n",
    "precip = pd.read_csv(\n",
    "    \"data/rain_nb1.csv\", parse_dates=[\"date\"], index_col=\"date\"\n",
    ").squeeze()\n",
    "print(\"The data type of the precip series is: %s\" % type(precip))\n",
    "\n",
    "# Import observed evaporation series\n",
    "evap = pd.read_csv(\n",
    "    \"data/evap_nb1.csv\", parse_dates=[\"date\"], index_col=\"date\"\n",
    ").squeeze()\n",
    "print(\"The data type of the evap series is: %s\" % type(evap))\n",
    "\n",
    "# Calculate the recharge to the groundwater\n",
    "recharge = precip - evap\n",
    "recharge.name = \"recharge\"  # set name if pandas series\n",
    "print(\"The data type of the recharge series is: %s\" % type(recharge))\n",
    "\n",
    "# Plot the time series of the precipitation and evaporation\n",
    "plt.figure()\n",
    "recharge.plot(label=\"Recharge\", figsize=(10, 4))\n",
    "plt.xlabel(\"Time [years]\")\n",
    "plt.ylabel(\"Recharge (m/year)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Create the time series model\n",
    "In this code block the actual time series model is created. First, an instance of the `Model` class is created (named `ml` here). Second, the different components of the time series model are created and added to the model. The imported time series are automatically checked for missing values and other inconsistencies. The keyword argument fillnan can be used to determine how missing values are handled. If any nan-values are found this will be reported by `pastas`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a model object by passing it the observed series\n",
    "ml = ps.Model(gwdata, name=\"GWL\")\n",
    "\n",
    "# Add the recharge data as explanatory variable\n",
    "sm = ps.StressModel(recharge, ps.Gamma(), name=\"recharge\", settings=\"evap\")\n",
    "ml.add_stressmodel(sm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Solve the model\n",
    "The next step is to compute the optimal model parameters. The default solver  uses a non-linear least squares method for the optimization. The python package `scipy` is used (info on `scipy's` least_squares solver can be found [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html)). Some standard optimization statistics are reported along with the optimized parameter values and correlations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Plot the results\n",
    "The solution can be plotted after a solution has been obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Advanced plotting\n",
    "There are many ways to further explore the time series model. `pastas` has some built-in functionalities that will provide the user with a quick overview of the model. The `plots` subpackage contains all the options. One of these is the method `plots.results` which provides a plot with more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.results(figsize=(10, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Statistics\n",
    "The `stats` subpackage includes a number of statistical functions that may applied to the model. One of them is the `summary` method, which gives a summary of the main statistics of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.stats.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Improvement: estimate evaporation factor\n",
    "In the previous model, the recharge was estimated as precipitation minus potential evaporation. A better model is to estimate the actual evaporation as a factor (called the evaporation factor here) times the potential evaporation. First, new model is created (called `ml2` here so that the original model `ml` does not get overwritten). Second, the `RechargeModel` object with a `Linear` recharge model is created, which combines the precipitation and evaporation series and adds a parameter for the evaporation factor `f`. The `RechargeModel` object is added to the model, the model is solved, and the results and statistics are plotted to the screen. Note that the new model gives a better fit (lower root mean squared error and higher explained variance), but that the Akiake information criterion indicates that the addition of the additional parameter does not improve the model signficantly (the Akaike criterion for model `ml2` is higher than for model `ml`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a model object by passing it the observed series\n",
    "ml2 = ps.Model(gwdata)\n",
    "\n",
    "# Add the recharge data as explanatory variable\n",
    "ts1 = ps.RechargeModel(\n",
    "    precip,\n",
    "    evap,\n",
    "    ps.Gamma(),\n",
    "    name=\"rainevap\",\n",
    "    recharge=ps.rch.Linear(),\n",
    "    settings=(\"prec\", \"evap\"),\n",
    ")\n",
    "ml2.add_stressmodel(ts1)\n",
    "\n",
    "# Solve the model\n",
    "ml2.solve()\n",
    "\n",
    "# Plot the results\n",
    "ml2.plot()\n",
    "\n",
    "# Statistics\n",
    "ml2.stats.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Origin of the series\n",
    "* The rainfall data is taken from rainfall station Heibloem in The Netherlands.\n",
    "* The evaporation data is taken from weather station Maastricht in The Netherlands.\n",
    "* The head data is well B58C0698, which was obtained from Dino loket"
   ]
  }
 ],
 "metadata": {
  "CodeCell": {
   "cm_config": {
    "lineWrapping": true
   }
  },
  "MarkdownCell": {
   "cm_config": {
    "lineWrapping": true
   }
  },
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "pastas_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.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "29475f8be425919747d373d827cb41e481e140756dd3c75aa328bf3399a0138e"
   }
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}