{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "999bc8b7-804c-4301-b8c7-5077a38dc06f",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Determining the crop coefficient function with Python\n",
    "\n",
    "*M. Vremec, October 2022, University of Graz*\n",
    "\n",
    "Data source: ZAMG - https://data.hub.zamg.ac.at\n",
    "\n",
    "What is done:\n",
    "\n",
    "- load the station data from ZAMG\n",
    "- estimate potential evapotranspiration\n",
    "- determine the crop coefficient function based on equation 65 in Allen et al. 1998\n",
    "- plot and store result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a2ae568-3c3b-44f5-8b82-ab5a99d35d2d",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pyet\n",
    "pyet.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da6f3e98-0073-4db9-b763-03beb2611e78",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 1. Loading daily data from ZAMG (Messstationen Tagesdaten)\n",
    "\n",
    "station: Graz Universität 16412\n",
    "\n",
    "Selected variables:\n",
    "- globalstrahlung (global radiation), J/cm2 needs to be in MJ/m3d, ZAMG abbreviation - strahl\n",
    "- arithmetische windgeschwindigkeit (wind speed), m/s, ZAMG abbreviation - vv\n",
    "- relative feuchte (relative humidity), %, ZAMG abbreviation - rel\n",
    "- lufttemparatur (air temperature) in 2 m, C, ZAMG abbreviation - t\n",
    "- lufttemperatur (air temperature) max in 2 m, C, ZAMG abbreviation - tmax\n",
    "- lufttemperatur (air temperature) min in 2 m, C, ZAMG abbreviation - tmin\n",
    "- latitute and elevation of a station"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "236434b2-3c33-4772-93ec-de0cb7317209",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "#read data\n",
    "data_16412 = pd.read_csv('data/example_1/klima_daily.csv', index_col=1, parse_dates=True)\n",
    "data_16412"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22ebf3bc-6b9a-4f81-a455-7ca1b51879dd",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 2. Calculate PET for Graz Universität - 16412"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc5ba933-22c2-4c5b-9ca6-43a4bcdad344",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Convert Glabalstrahlung J/cm2 to MJ/m2 by dividing to 100\n",
    "\n",
    "meteo = pd.DataFrame({\"time\":data_16412.index, \"tmean\":data_16412.t, \"tmax\":data_16412.tmax, \"tmin\":data_16412.tmin, \"rh\":data_16412.rel, \n",
    "                      \"wind\":data_16412.vv, \"rs\":data_16412.strahl/100})\n",
    "time, tmean, tmax, tmin, rh, wind, rs = [meteo[col] for col in meteo.columns]\n",
    "\n",
    "lat = 47.077778*np.pi/180  # Latitude of the meteorological station, converting from degrees to radians\n",
    "elevation = 367  # meters above sea-level\n",
    "\n",
    "# Estimate potential ET with Penman-Monteith FAO-56\n",
    "pet_pm = pyet.pm_fao56(tmean, wind, rs=rs, elevation=elevation, \n",
    "                      lat=lat, tmax=tmax, tmin=tmin, rh=rh)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "696cb6aa-773d-4e1b-8a11-64ee5a5242e9",
   "metadata": {},
   "source": [
    "## 3. Determine the crop coefficient function\n",
    "\n",
    "Based on: https://www.fao.org/3/x0490e/x0490e0b.htm\n",
    "figure 34.\n",
    "\n",
    "\n",
    "![Figure 34](https://www.fao.org/3/x0490e/x0490e6k.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92b3767-7233-461c-a924-cafc6f5c9882",
   "metadata": {},
   "outputs": [],
   "source": [
    "Kcini = 0.3 \n",
    "Kcmid = 1.1\n",
    "Kcend = 0.65\n",
    "\n",
    "crop_ini = pd.Timestamp(\"2020-04-01\")\n",
    "crop_dev = pd.Timestamp(\"2020-05-01\")\n",
    "mid_season = pd.Timestamp(\"2020-06-01\")\n",
    "late_s_start = pd.Timestamp(\"2020-07-01\")\n",
    "late_s_end = pd.Timestamp(\"2020-08-01\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86d39fcb-0dea-42ca-9b7a-afa5327ad990",
   "metadata": {},
   "outputs": [],
   "source": [
    "kc = pd.Series(index=[crop_ini, crop_dev, mid_season, late_s_start, late_s_end],\n",
    "                 data=[Kcini, Kcini, Kcmid, Kcmid, Kcend])\n",
    "kc = kc.resample(\"d\").mean().interpolate()\n",
    "kc.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c73067c3-2d7e-4bb4-8cff-91f707f7c34b",
   "metadata": {},
   "outputs": [],
   "source": [
    "petc = pet_pm.loc[crop_dev:late_s_end] * kc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e2564c-48ed-46b7-a275-fd17db592456",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## 4. Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9a4582f-219b-47f8-9ed2-dac80720cf55",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "pet_pm.loc[crop_dev:late_s_end].plot(label=\"Potential evapotranspiration\")\n",
    "petc.loc[crop_dev:late_s_end].plot(label=\"Potential crop evapotranspiration\")"
   ]
  }
 ],
 "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": 5
}