{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "999bc8b7-804c-4301-b8c7-5077a38dc06f",
   "metadata": {},
   "source": [
    "# Potential Evapotranspiration from ZAMG INCA data (NetCDF)\n",
    "*M. Vremec, October 2022, University of Graz*\n",
    "\n",
    "\n",
    "What is done:\n",
    "\n",
    "- load the data from ZAMG\n",
    "- estimate potential evapotranspiration\n",
    "- plot and store results\n",
    "\n",
    "Data source: ZAMG - https://data.hub.zamg.ac.at"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a2ae568-3c3b-44f5-8b82-ab5a99d35d2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "#import netcdf4  # Needs to be installed\n",
    "import pyet\n",
    "pyet.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da6f3e98-0073-4db9-b763-03beb2611e78",
   "metadata": {},
   "source": [
    "## 1. Loading daily data from ZAMG (INCA hourly NetCDF data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "236434b2-3c33-4772-93ec-de0cb7317209",
   "metadata": {},
   "outputs": [],
   "source": [
    "#read data\n",
    "xr_ds = xr.open_dataset(\"data/example_2/incal_hourly_20120501T0000_20120930T2300.nc\", \n",
    "            engine=\"netcdf4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd56788b-5aed-499b-8a6a-9c6fb68bb2f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resample and define input meteorological variables\n",
    "tmean = xr_ds[\"T2M\"].resample(time=\"1D\").mean()\n",
    "tmax = xr_ds[\"T2M\"].resample(time=\"1D\").max()\n",
    "tmin = xr_ds[\"T2M\"].resample(time=\"1D\").min()\n",
    "rh = xr_ds[\"RH2M\"].resample(time=\"1D\").mean()\n",
    "rhmax = xr_ds[\"RH2M\"].resample(time=\"1D\").max()\n",
    "rhmin = xr_ds[\"RH2M\"].resample(time=\"1D\").min()\n",
    "wind = ((np.abs(xr_ds[\"VV\"]) + np.abs(xr_ds[\"UU\"])) / 2).resample(time=\"1D\").mean()\n",
    "rs = xr_ds[\"GL\"].resample(time=\"1D\").mean() * 86400 / 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5bb9668-075b-4a8c-99e5-64a998c9c9d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define latitude and elevation\n",
    "lat =  tmean.lat * np.pi / 180  \n",
    "elevation = lat / lat * 350"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22ebf3bc-6b9a-4f81-a455-7ca1b51879dd",
   "metadata": {},
   "source": [
    "## 2. Calculate PET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "542a2a7a-7354-4ce7-a7d7-644f8f650238",
   "metadata": {},
   "outputs": [],
   "source": [
    "lat1 = 80 * np.pi / 180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b8859a2-37e3-40db-afc1-01b449adfc35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate evapotranspiration with nine different methods \n",
    "pet_penman = pyet.penman(tmean, wind, rs=rs, elevation=elevation, lat=lat, tmax=tmax, tmin=tmin, rh=rh)\n",
    "pet_pt = pyet.priestley_taylor(tmean, rs=rs, elevation=elevation, lat=lat, tmax=tmax, tmin=tmin, rh=rh)\n",
    "pet_makkink = pyet.makkink(tmean, rs, elevation=elevation)\n",
    "pet_fao56 = pyet.pm_fao56(tmean, wind, rs=rs, elevation=elevation, lat=lat, tmax=tmax, tmin=tmin, rh=rh)\n",
    "pet_hamon = pyet.hamon(tmean, lat=lat, method=1)\n",
    "pet_oudin = pyet.oudin(tmean, lat=lat)\n",
    "pet_haude = pyet.haude(tmax, rh)\n",
    "pet_turc = pyet.turc(tmean, rs, rh)\n",
    "pet_har = pyet.hargreaves(tmean, tmax, tmin, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e2564c-48ed-46b7-a275-fd17db592456",
   "metadata": {},
   "source": [
    "## 3. Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c613f4-edc1-4472-8ab3-5ed0e1d373f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(ncols=2, figsize=(12,3))\n",
    "pet_penman[:,2,2].plot(ax=axs[0], label=\"Penman\")\n",
    "pet_pt[:,2,2].plot(ax=axs[0], label=\"Priestley-Taylor\")\n",
    "pet_makkink[:,2,2].plot(ax=axs[0], label=\"Makkink\")\n",
    "pet_fao56[:,2,2].plot(ax=axs[0], label=\"FAO56\")\n",
    "pet_hamon[:,2,2].plot(ax=axs[0], label=\"Hamon\")\n",
    "pet_oudin[:,2,2].plot(ax=axs[0], label=\"Oudin\")\n",
    "pet_haude[:,2,2].plot(ax=axs[0], label=\"Haude\")\n",
    "pet_turc[:,2,2].plot(ax=axs[0], label=\"Turc\")\n",
    "pet_har[:,2,2].plot(ax=axs[0], label=\"Hargreaves\")\n",
    "axs[0].legend()\n",
    "\n",
    "pet_penman[:,2,2].cumsum().plot(ax=axs[1], label=\"Penman\")\n",
    "pet_pt[:,2,2].cumsum().plot(ax=axs[1], label=\"Priestley-Taylor\")\n",
    "pet_makkink[:,2,2].cumsum().plot(ax=axs[1], label=\"Makkink\")\n",
    "pet_fao56[:,2,2].cumsum().plot(ax=axs[1], label=\"FAO56\")\n",
    "pet_hamon[:,2,2].cumsum().plot(ax=axs[1], label=\"Hamon\")\n",
    "pet_oudin[:,2,2].cumsum().plot(ax=axs[1], label=\"Oudin\")\n",
    "pet_haude[:,2,2].cumsum().plot(ax=axs[1], label=\"Haude\")\n",
    "pet_turc[:,2,2].cumsum().plot(ax=axs[1], label=\"Turc\")\n",
    "pet_har[:,2,2].cumsum().plot(ax=axs[1], label=\"Hargreaves\")\n",
    "axs[1].legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6697f56-83ef-49f4-8037-9235f6e88a51",
   "metadata": {},
   "source": [
    "## 4. Compare point with pandas.Series vs. point from xarray.DataArray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b70386d-18f3-4485-bdc1-8969630018b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmeans = xr_ds[\"T2M\"].resample(time=\"1D\").mean()[:, 4, 4].to_series()\n",
    "tmaxs = xr_ds[\"T2M\"].resample(time=\"1D\").max()[:, 4, 4].to_series()\n",
    "tmins = xr_ds[\"T2M\"].resample(time=\"1D\").min()[:, 4, 4].to_series()\n",
    "rhs = xr_ds[\"RH2M\"].resample(time=\"1D\").mean()[:, 4, 4].to_series()\n",
    "rhmaxs = xr_ds[\"RH2M\"].resample(time=\"1D\").max()[:, 4, 4].to_series()\n",
    "rhmins = xr_ds[\"RH2M\"].resample(time=\"1D\").min()[:, 4, 4].to_series()\n",
    "winds = ((np.abs(xr_ds[\"VV\"]) + np.abs(xr_ds[\"UU\"])) / 2).resample(time=\"1D\").mean()[:, 4, 4].to_series()\n",
    "rss = (xr_ds[\"GL\"].resample(time=\"1D\").mean() * 86400 / 1000000)[:, 4, 4].to_series()\n",
    "lats = float(lat[4, 4])\n",
    "elevations = float(elevation[4, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef26699-5611-4c93-bb26-4fdef1138eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate evapotranspiration with nine different methods with Pandas.Series\n",
    "pet_penmans = pyet.penman(tmeans, winds, rs=rss, elevation=elevations, lat=lats, tmax=tmaxs, tmin=tmins, rh=rhs)\n",
    "pet_pts = pyet.priestley_taylor(tmeans, rs=rss, elevation=elevations, lat=lats, tmax=tmaxs, tmin=tmins, rh=rhs)\n",
    "pet_makkinks = pyet.makkink(tmeans, rss, elevation=elevations)\n",
    "pet_fao56s = pyet.pm_fao56(tmeans, winds, rs=rss, elevation=elevations, lat=lats, tmax=tmaxs, tmin=tmins, rh=rhs)\n",
    "pet_hamons = pyet.hamon(tmeans, lat=lats, method=1)\n",
    "pet_oudins = pyet.oudin(tmeans, lat=lats)\n",
    "pet_haudes = pyet.haude(tmaxs, rhs)\n",
    "pet_turcs = pyet.turc(tmeans, rss, rhs)\n",
    "pet_hars = pyet.hargreaves(tmeans, tmaxs, tmins, lats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7165b896-c86e-460f-a592-874f2381f7d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fm = \"\"\"\n",
    "    ABC\n",
    "    DEF\n",
    "    GHI\n",
    "    \"\"\"\n",
    "fig,axs = plt.subplot_mosaic(mosaic=fm, figsize=(12,8))\n",
    "axs[\"A\"].scatter(pet_penman[:,4,4].values, pet_penmans.values)\n",
    "axs[\"B\"].scatter(pet_pt[:,4,4].values, pet_pts.values)\n",
    "axs[\"C\"].scatter(pet_makkink[:,4,4].values, pet_makkinks.values)\n",
    "axs[\"D\"].scatter(pet_fao56[:,4,4].values, pet_fao56s.values)\n",
    "axs[\"E\"].scatter(pet_hamon[:,4,4].values, pet_hamons.values)\n",
    "axs[\"F\"].scatter(pet_oudin[:,4,4].values, pet_oudins.values)\n",
    "axs[\"G\"].scatter(pet_haude[:,4,4].values, pet_haudes.values)\n",
    "axs[\"H\"].scatter(pet_turc[:,4,4].values, pet_turcs.values)\n",
    "axs[\"I\"].scatter(pet_har[:,4,4].values, pet_hars.values)\n",
    "for i in axs.keys():\n",
    "    axs[i].plot([0,7], [0,7])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fe2a55b-838a-4e3a-a621-c5fe55db399a",
   "metadata": {},
   "source": [
    "## 5. Store results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c5e47fb-3ec7-4bc9-ad6c-be0edc67ba09",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pet_pt.to_netcdf('../pe_pt_INCA.csv')"
   ]
  }
 ],
 "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
}