{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibration of potential evapotranspiration methods\n",
    "*M. Vremec, R.A. Collenteur University of Graz, 2021*\n",
    "\n",
    "In this notebook it is shown how to calibrate various (potential) evapotranspiration (PET) equations, using a linear regression relationship between daily potential evapotranspiration obtained with the FAO-56 equation and daily PET estimates obtained with the alternative methods. This notebook requires Scipy and SkLearn python packages to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import least_squares\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "import pyet as pyet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load KNMI Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"data/example_3/etmgeg_260.txt\", skiprows=46, delimiter=\",\", \n",
    "                   skipinitialspace=True, index_col=\"YYYYMMDD\", parse_dates=True).loc[\"2018\",:]\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.Estimation of potential evapotranspiration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define meteorological variables needed for PE estimation\n",
    "meteo = pd.DataFrame({\"tmean\":data.TG/10, \"tmax\":data.TX/10, \"tmin\":data.TN/10, \"rh\":data.UG, \"wind\":data.FG/10, \"rs\":data.Q/100})\n",
    "tmean, tmax, tmin, rh, wind, rs = [meteo[col] for col in meteo.columns]\n",
    "pressure =  data.PG / 100 # to kPa\n",
    "wind = data.FG/10  # to m/s\n",
    "lat = 0.91  # [rad]\n",
    "elevation = 4 \n",
    "\n",
    "pet_fao56 = pyet.pm_fao56(tmean, wind, rs=rs, elevation=elevation, lat=lat, \n",
    "                         tmax=tmax, tmin=tmin, rh=rh)  # FAO-56 method\n",
    "pet_romanenko = pyet.romanenko(tmean, rh)  # Romanenko method\n",
    "pet_abtew = pyet.abtew(tmean, rs)  # Abtew method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model performance of the Abtew and Romanenko method will be evaluated using the root mean square error (RMSE), as implemented in SkLearn’s mean_squared_error method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pet_fao56.plot()\n",
    "pet_romanenko.plot()\n",
    "pet_abtew.plot()\n",
    "plt.ylabel(\"PE [mm/day]\")\n",
    "print(\"RMSE(Romanenko) = {} mm/d\".format(mean_squared_error(pet_fao56, pet_romanenko, squared=False)))\n",
    "print(\"RMSE(Abtew) = {} mm/d\".format(mean_squared_error(pet_fao56, pet_abtew, squared=False)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Calibration of different PE equations\n",
    "\n",
    "The least squares approach is applied to estimate the parameters in the PE equations, by minimizing the sum of the residuals between the simulated (Abtew and Romanenko) and observed (FAO-56) data. The minimization of the objective function is done using the Trust Region Reflective algorithm, as implemented in Scipy’s least squares method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Romanenko method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function for computing residuals\n",
    "def cal_romanenko(k, obs):\n",
    "    return pyet.romanenko(tmean, rh, k)-obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimate k in the Romanenko method\n",
    "x0 = 4.5  # initial estimate of parameter\n",
    "res_1 = least_squares(cal_romanenko, x0, args=[pet_fao56])\n",
    "res_1.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute RMSE using the calibrated value of k\n",
    "pet_romanenko_cal = pyet.romanenko(tmean, rh, k=res_1.x)\n",
    "print(\"RMSE(Romanenko) = {} mm/d\".format(mean_squared_error(pet_fao56, pet_romanenko_cal, squared=False)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RMSE (calibrated) = 0.546 < RMSE (uncalibrated) = 0.694"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Abtew method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function for computing residuals and initial estimate of parameters\n",
    "def cal_abtew(k,obs):\n",
    "    return pyet.abtew(tmean, rs, k)-obs\n",
    "x0 = 0.53"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimate k in the Romanenko method\n",
    "res_2 = least_squares(cal_abtew, x0, args=[pet_fao56])\n",
    "res_2.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pet_abtew_cali = pyet.abtew(tmean, rs, res_2.x)\n",
    "print(\"RMSE(Abtew) = {} mm/d\".format(mean_squared_error(pet_fao56, pet_abtew_cali, squared=False)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RMSE (calibrated) = 0.613 < RMSE (uncalibrated) = 0.741"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pet_fao56.plot()\n",
    "pet_romanenko_cal.plot()\n",
    "pet_abtew_cali.plot()\n",
    "plt.ylabel(\"PET [mm/d]\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}