{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d0959a4c",
   "metadata": {},
   "source": [
    "# Potential Evapotranspiration from KNMI data\n",
    "*R.A. Collenteur, Eawag, 2023*\n",
    "\n",
    "Data source: KNMI - https://dataplatform.knmi.nl/\n",
    "\n",
    "In this notebook it is shown how to compute (potential) evapotranspiration from meteorological data using PyEt. Meteorological data is observed by the KNMI at De Bilt in the Netherlands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8c6cab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pyet\n",
    "pyet.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61cb9b9e",
   "metadata": {},
   "source": [
    "## 1. Load KNMI Data\n",
    "\n",
    "We first load the raw meteorological data observed by the KNMI at De Bilt in the Netherlands. This datafile contains a lot of different variables, please see the end of the notebook for an explanation of all the variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbc6c337",
   "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",
   "id": "be888ed3",
   "metadata": {},
   "source": [
    "## 2. Estimating potential evapotranspiration\n",
    "\n",
    "Now that we have the input data, we can estimate potential evapotranspiration with different estimation methods. Here we choose the Penman, Priestley-Taylor, Makkink, and Oudin methods. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abe8a230",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preprocess the input data\n",
    "meteo = pd.DataFrame({\"tmean\":data.TG/10, \"tmax\":data.TX/10, \"tmin\":data.TN/10, \n",
    "                      \"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  # Latitude of the meteorological station\n",
    "elevation = 4  # meters above sea-level \n",
    "\n",
    "# Estimate evapotranspiration with four different methods\n",
    "pet_penman = pyet.penman(tmean, wind, rs=rs, elevation=4, lat=0.91, tmax=tmax, tmin=tmin, rh=rh)\n",
    "pet_pt = pyet.priestley_taylor(tmean, rs=rs, elevation=4, lat=0.91, tmax=tmax, tmin=tmin, rh=rh)\n",
    "pet_makkink = pyet.makkink(tmean, rs, elevation=4, pressure=pressure)\n",
    "pet_oudin = pyet.oudin(tmean, lat=0.91)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578f8ac2",
   "metadata": {},
   "source": [
    "## 3. Plot the results\n",
    "\n",
    "We plot the cumulative sums to compare the different estimation methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c5a8008",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(figsize=(14,3.5), ncols=2)\n",
    "\n",
    "pet = [pet_penman, pet_pt, pet_makkink, pet_makkink]\n",
    "names = [\"Penman\", \"Priestley-Taylor\", \"Makkink - PyEt\", \"Oudin\"]\n",
    "\n",
    "for df, name in zip(pet, names):\n",
    "    axs[0].plot(df,label=name)\n",
    "    axs[1].plot(df.cumsum(),label=name)\n",
    "\n",
    "axs[0].set_ylabel(\"PET [mm/day]\", fontsize=16)\n",
    "axs[1].set_ylabel(\"Cumulative PET [mm]\", fontsize=12)\n",
    "\n",
    "for i in (0,1):\n",
    "    axs[i].set_xlabel(\"Date\", fontsize=12)\n",
    "    axs[i].legend(loc=2, fontsize=14)\n",
    "    axs[i].tick_params(\"both\", direction=\"in\", labelsize=14)\n",
    "plt.tight_layout()\n",
    "\n",
    "#plt.savefig(\"Figure1.png\", dpi=300)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "75dacdc4",
   "metadata": {},
   "source": [
    "## 4. Comparison: pyet Makkink vs KNMI Makkink\n",
    "\n",
    "The KNMI also provides Makkink potential evapotranspiration data (column EV24). We can now compare the results from Pyet and the KNMI to confirm that these are roughly the same. Pyet also has a method `makkink_knmi` that calculates exactly the same Makkink evaporation as the KNMI but that is only suitable for the Netherlands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95c6b234",
   "metadata": {},
   "outputs": [],
   "source": [
    "pet_knmi = data.EV24 / 10  # Makkink potential evaporation computen by the KNMI for comparison\n",
    "pet_makkink_knmi = pyet.makkink_knmi(tmean, rs).round(1) # same as pet_knmi (if rounded up to 1 decimal) but calculated from tmean and rs\n",
    "\n",
    "# Plot the two series against each other\n",
    "_, ax = plt.subplots(1, 2, figsize=(9,4), sharex=True, sharey=True)\n",
    "ax[0].scatter(pet_makkink, pet_knmi, s=15, color=\"C0\")\n",
    "ax[0].plot([0,6],[0,6], color=\"red\", label=\"1:1 line\")\n",
    "ax[0].set_xlabel(\"pyet-Makkink 'EV24' [mm]\")\n",
    "ax[0].grid()\n",
    "\n",
    "ax[1].scatter(pet_makkink_knmi, pet_knmi, s=15, color=\"C1\")\n",
    "ax[1].plot([0,6],[0,6], color=\"red\")\n",
    "ax[1].set_xlabel(\"pyet-Makkink_KNMI [mm]\")\n",
    "ax[1].grid()\n",
    "\n",
    "ax[0].set_ylabel(\"KNMI-Makkink [mm]\")\n",
    "ax[0].legend(loc=(0,1), frameon=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51c36d1",
   "metadata": {},
   "source": [
    "## 5. Estimation of PET - all methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a2b17d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "pet_df = pyet.calculate_all(tmean, wind, rs, elevation, lat, tmax=tmax,\n",
    "                                tmin=tmin, rh=rh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "055e90d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(figsize=(13,4), ncols=2)\n",
    "pet_df.plot(ax=axs[0])\n",
    "pet_df.cumsum().plot(ax=axs[1], legend=False)\n",
    "\n",
    "axs[0].set_ylabel(\"PET [mm/day]\", fontsize=12)\n",
    "axs[1].set_ylabel(\"Cumulative PET [mm]\", fontsize=12)\n",
    "axs[0].legend(ncol=6, loc=[0,1.])\n",
    "for i in (0,1):\n",
    "    axs[i].set_xlabel(\"Date\", fontsize=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e55de8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.savefig(\"PET_methods.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "815addc4",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Column description of the KNMI data\n",
    "\n",
    "- DDVEC     = Vector mean wind direction in degrees (360=north, 90=east, 180=south, 270=west, 0=calm/variable)\n",
    "- FHVEC     = Vector mean windspeed (in 0.1 m/s)\n",
    "- FG        = Daily mean windspeed (in 0.1 m/s) \n",
    "- FHX       = Maximum hourly mean windspeed (in 0.1 m/s)\n",
    "- FHXH      = Hourly division in which FHX was measured\n",
    "- FHN       = Minimum hourly mean windspeed (in 0.1 m/s)\n",
    "- FHNH      = Hourly division in which FHN was measured\n",
    "- FXX       = Maximum wind gust (in 0.1 m/s)\n",
    "- FXXH      = Hourly division in which FXX was measured\n",
    "- TG        = Daily mean temperature in (0.1 degrees Celsius)\n",
    "- TN        = Minimum temperature (in 0.1 degrees Celsius)\n",
    "- TNH       = Hourly division in which TN was measured\n",
    "- TX        = Maximum temperature (in 0.1 degrees Celsius)\n",
    "- TXH       = Hourly division in which TX was measured\n",
    "- T10N      = Minimum temperature at 10 cm above surface (in 0.1 degrees Celsius)\n",
    "- T10NH     = 6-hourly division in which T10N was measured; 6=0-6 UT, 12=6-12 UT, 18=12-18 UT, 24=18-24 UT \n",
    "- SQ        = Sunshine duration (in 0.1 hour) calculated from global radiation (-1 for <0.05 hour)\n",
    "- SP        = Percentage of maximum potential sunshine duration\n",
    "- Q         = Global radiation (in J/cm2)\n",
    "- DR        = Precipitation duration (in 0.1 hour)\n",
    "- RH        = Daily precipitation amount (in 0.1 mm) (-1 for <0.05 mm)\n",
    "- RHX       = Maximum hourly precipitation amount (in 0.1 mm) (-1 for <0.05 mm)\n",
    "- RHXH      = Hourly division in which RHX was measured\n",
    "- PG        = Daily mean sea level pressure (in 0.1 hPa) calculated from 24 hourly values\n",
    "- PX        = Maximum hourly sea level pressure (in 0.1 hPa)\n",
    "- PXH       = Hourly division in which PX was measured\n",
    "- PN        = Minimum hourly sea level pressure (in 0.1 hPa)\n",
    "- PNH       = Hourly division in which PN was measured\n",
    "- VVN       = Minimum visibility; 0: <100 m, 1:100-200 m, 2:200-300 m,..., 49:4900-5000 m, 50:5-6 km, 56:6-7 km, 57:7-8 km,..., 79:29-30 km, 80:30-35 km, 81:35-40 km,..., 89: >70 km)\n",
    "- VVNH      = Hourly division in which VVN was measured\n",
    "- VVX       = Maximum visibility; 0: <100 m, 1:100-200 m, 2:200-300 m,..., 49:4900-5000 m, 50:5-6 km, 56:6-7 km, 57:7-8 km,..., 79:29-30 km, 80:30-35 km, 81:35-40 km,..., 89: >70 km)\n",
    "- VVXH      = Hourly division in which VVX was measured\n",
    "- NG        = Mean daily cloud cover (in octants, 9=sky invisible)\n",
    "- UG        = Daily mean relative atmospheric humidity (in percents)\n",
    "- UX        = Maximum relative atmospheric humidity (in percents)\n",
    "- UXH       = Hourly division in which UX was measured\n",
    "- UN        = Minimum relative atmospheric humidity (in percents)\n",
    "- UNH       = Hourly division in which UN was measured\n",
    "- EV24      = Potential evapotranspiration (Makkink) (in 0.1 mm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df028060-9d5f-4616-9797-d44e6b66002b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  },
  "vscode": {
   "interpreter": {
    "hash": "8a365f246fd2b6cc3f433897ea78f6113fbca3a5aebe59945b06a3e6c75c5dde"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}