{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read KNMI observations using hydropandas\n",
    "\n",
    "This notebook introduces how to use the `hydropandas` package to read, process and visualise KNMI data.\n",
    "\n",
    "Martin & Onno - 2022"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=top></a>Notebook contents\n",
    "\n",
    "1. [Observation types](#Obs)\n",
    "2. [Get KNMI data](#reading)\n",
    "3. [Get ObsCollections](#readingOC)\n",
    "4. [Precipitation data](#precipitation)\n",
    "5. [Reference evaporation types](#EvapRef)\n",
    "6. [Spatial interpolation of evaporation](#EvapInterp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from IPython.display import display\n",
    "from scipy.interpolate import NearestNDInterpolator, RBFInterpolator\n",
    "from tqdm import tqdm\n",
    "\n",
    "import hydropandas as hpd\n",
    "from hydropandas.io import knmi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hpd.util.get_color_logger(\"INFO\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Observation types](#top)<a id=Obs></a>\n",
    "\n",
    "The hydropandas package has a function to read all kinds of KNMI observations. These are stored in an `Obs` object. There are three types of observations you can obtain from the KNMI:\n",
    "\n",
    "- `EvaporationObs`, for evaporation time series\n",
    "- `PrecipitationObs`, for precipitation time series\n",
    "- `MeteoObs`, for all the other meteorological time series\n",
    "\n",
    "Get evaporation in [m/day] for KNMI station 344 (Rotterdam Airport)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o = hpd.EvaporationObs.from_knmi(stn=344, start=\"2022\")\n",
    "display(o.head())\n",
    "o.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get precipitation in [m/day] for KNMI station 344 (Rotterdam Airport)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o = hpd.PrecipitationObs.from_knmi(stn=344, start=\"2022\")\n",
    "display(o.head())\n",
    "o.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get temperature in [°C] for KNMI station 344 (Rotterdam Airport)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o = hpd.MeteoObs.from_knmi(stn=344, meteo_var=\"TG\", start=\"2022\")\n",
    "display(o.head())\n",
    "o.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**attributes**\n",
    "\n",
    "All `PrecipitationObs`, `EvaporationObs` and `MeteoObs` objects have the attributes:\n",
    "\n",
    "* `name`: meteo variable and station name\n",
    "* `x`: x-coordinate in m RD\n",
    "* `y`: y-coordinate in m RD\n",
    "* `station`: station number\n",
    "* `unit`: measurement unit\n",
    "* `meta`: dictionary with other metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"name: {o.name}\")\n",
    "print(f\"x,y: {(o.x, o.y)}\")\n",
    "print(f\"station: {o.station}\")\n",
    "print(\"unit\", o.unit)\n",
    "print(\"metadata:\")\n",
    "for key, item in o.meta.items():\n",
    "    print(f\"    {key}: {item}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Get KNMI data](#top)<a id=reading></a>\n",
    "\n",
    "The `from_knmi` method can be used to obtain one type of measurements at one station. There are two ways to select the station.\n",
    "\n",
    "1. with the station number `stn`\n",
    "2. by specifying the `xy` coordinates, the station nearest to this coordinates is used.\n",
    "\n",
    "Below both methods are used to obtain precipitation measurements. Notice that they return the same data because station 344 is nearest to the given xy coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o1 = hpd.PrecipitationObs.from_knmi(stn=344)\n",
    "o2 = hpd.PrecipitationObs.from_knmi(xy=(90600, 442800))\n",
    "o1.equals(o2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**read options**\n",
    "\n",
    "The `from_knmi` method has other optional arguments:\n",
    "\n",
    "- `stn`: station number.\n",
    "- `start`: the start date of the time series you want, default is 1st of January 2019.\n",
    "- `end`: the end date of the time series you want, default is today.\n",
    "- `fill_missing_obs`: option to fill missing values with values from the nearest KNMI station. If measurements are filled an extra column is added to the time series in which the station number is shown that was used to fill a particular missing value.\n",
    "- `interval`: time interval of the time series, default is 'daily'\n",
    "- `use_api`: select the method to obtain the data, default is 'True'\n",
    "- `raise_exception`: option to raise an error when the requested time series is empty.\n",
    "***\n",
    "\n",
    "The 3 examples below give a brief summary of these options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example 1 get daily average temperature (TG) from 1900 till now\n",
    "o_t = hpd.MeteoObs.from_knmi(\"TG\", stn=344, start=\"1960\")\n",
    "o_t.plot(figsize=(16, 4), grid=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example 2 get daily average precipitation from 1972 with and without filling missing measurements\n",
    "f, axes = plt.subplots(figsize=(16, 6), nrows=2, sharex=True, sharey=True)\n",
    "\n",
    "o_rd = hpd.PrecipitationObs.from_knmi(\n",
    "    \"RD\", stn=892, start=\"1972\", end=\"2023\", fill_missing_obs=False\n",
    ")\n",
    "o_rd.plot(figsize=(16, 4), ax=axes[0])\n",
    "\n",
    "o_rd_filled = hpd.PrecipitationObs.from_knmi(\n",
    "    \"RD\", stn=892, start=\"1972\", end=\"2023\", fill_missing_obs=True\n",
    ")\n",
    "\n",
    "o_rd_filled.loc[o_rd_filled[\"station\"] == \"892\", \"RD\"].plot(\n",
    "    ax=axes[1], label=\"oorspronkelijke metingen\"\n",
    ")\n",
    "o_rd_filled.loc[o_rd_filled[\"station\"] != \"892\", \"RD\"].plot(\n",
    "    ax=axes[1], color=\"orange\", label=\"opgevulde metingen\"\n",
    ")\n",
    "axes[0].set_title(\"precipitation Giersbergen (1972-2023)\")\n",
    "axes[1].legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see the station_opvulwaarde\n",
    "display(o_rd_filled.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example 3 get evaporation\n",
    "logging.getLogger().getEffectiveLevel()\n",
    "logging.getLogger().setLevel(logging.INFO)\n",
    "\n",
    "o_ev = hpd.EvaporationObs.from_knmi(\n",
    "    stn=344, start=\"1972\", end=\"2023\", fill_missing_obs=True\n",
    ")\n",
    "o_ev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Get ObsCollections](#top)<a id=readingOC></a>\n",
    "\n",
    "You can read multiple `Observation` objects at once using the `hpd.read_knmi` method. The observations are stored in an `ObsCollection` object. Below an example to obtain precipitation (RH) and evaporation (EV24) from the measurement stations Rotterdam (344) and De Bilt (260). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oc = hpd.read_knmi(stns=[344, 260], meteo_vars=[\"RH\", \"EV24\"])\n",
    "oc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "other options in the `read_knmi` function are:\n",
    "\n",
    "- specify `locations` as a dataframe with x, y coordinates (RD_new), the function will find the nearest KNMI station for every location.\n",
    "- specify `xmid` and `ymid` which are 2 arrays corresponding to a structured grid to obtain the nearest KNMI station for every cell in the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations = pd.DataFrame(index=[\"Rotterdam\"], data={\"x\": 90600, \"y\": 442800})\n",
    "display(locations)\n",
    "hpd.read_knmi(locations=locations, meteo_vars=[\"RH\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hpd.read_knmi(xy=((90600, 442800),), meteo_vars=[\"RH\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Precipitation](#top)<a id=precipitation></a>\n",
    "\n",
    "Using hydropandas 3 different precipitation products can be accessed:\n",
    "\n",
    "1. Daily data from a meteorological station\n",
    "2. Daily data from a neerslag (precipitation) station\n",
    "3. Hourly data from a meteorological station\n",
    "\n",
    "All three products can be obtained using the `from_knmi` method. Product 1 and 2 can also be accessed without the api.\n",
    "\n",
    "By default data is used from a meteorological station. Specify `meteo_var='RD'` to access data from a neerslag (precipitation) station."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = \"2010-1-1\"\n",
    "end = \"2010-1-10\"\n",
    "\n",
    "# daily meteo station (default)\n",
    "precip1 = hpd.PrecipitationObs.from_knmi(stn=260, start=start, end=end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# daily neerslag station\n",
    "precip2 = hpd.PrecipitationObs.from_knmi(meteo_var=\"RD\", stn=550, start=start, end=end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hourly meteo station (only works with api)\n",
    "precip3 = hpd.PrecipitationObs.from_knmi(\n",
    "    stn=260,\n",
    "    start=start,\n",
    "    end=end,\n",
    "    interval=\"hourly\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# daily meteo station without api\n",
    "precip4 = hpd.PrecipitationObs.from_knmi(\n",
    "    stn=260,\n",
    "    start=start,\n",
    "    end=end,\n",
    "    use_api=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# daily neerslag station without api\n",
    "precip5 = hpd.PrecipitationObs.from_knmi(\n",
    "    meteo_var=\"RD\",\n",
    "    stn=550,\n",
    "    start=start,\n",
    "    end=end,\n",
    "    use_api=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are the differences between the precipitation estimates from different station types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "precip1[\"RH\"].plot(\n",
    "    ax=ax,\n",
    "    drawstyle=\"steps\",\n",
    "    ls=\"--\",\n",
    "    lw=3,\n",
    "    label=str(precip1.station) + \"_dagelijks\",\n",
    ")\n",
    "\n",
    "precip2[\"RD\"].plot(\n",
    "    ax=ax,\n",
    "    drawstyle=\"steps\",\n",
    "    ls=\"--\",\n",
    "    lw=3,\n",
    "    label=str(precip2.station) + \"_dagelijks\",\n",
    ")\n",
    "\n",
    "precip3[\"RH\"].plot(\n",
    "    ax=ax,\n",
    "    drawstyle=\"steps\",\n",
    "    label=str(precip3.station) + \"_uurlijks\",\n",
    ")\n",
    "\n",
    "precip4[\"RH\"].plot(\n",
    "    ax=ax,\n",
    "    drawstyle=\"steps\",\n",
    "    marker=\"o\",\n",
    "    label=str(precip4.station) + \"_dagelijks_geen_api\",\n",
    ")\n",
    "\n",
    "precip5[\"RD\"].plot(\n",
    "    ax=ax,\n",
    "    drawstyle=\"steps\",\n",
    "    marker=\"o\",\n",
    "    label=str(precip5.station) + \"_dagelijks_geen_api\",\n",
    ")\n",
    "\n",
    "\n",
    "ax.legend()\n",
    "ax.grid()\n",
    "ax.set_ylabel(\"neerslag [m/dag]\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The locations of the stations can be plotted onto a map using the contextily package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import contextily as cx\n",
    "\n",
    "oc = hpd.ObsCollection([precip1, precip2])\n",
    "\n",
    "gdf = oc.to_gdf()\n",
    "\n",
    "gdf = gdf.set_crs(28992)\n",
    "\n",
    "gdf[\"name\"] = gdf.index\n",
    "\n",
    "ax = gdf.buffer(2000).plot(alpha=0, figsize=(8, 8))\n",
    "\n",
    "gdf.plot(\"name\", ax=ax, cmap=\"jet\", legend=True, markersize=100)\n",
    "\n",
    "\n",
    "cx.add_basemap(ax, crs=28992, source=cx.providers.OpenStreetMap.Mapnik)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Evaporation](#top)<a id=EvapRef></a>\n",
    "\n",
    "KNMI provides the Makking reference evaporation (meteo_var EV24). Hydropandas provides a posssibility to calculate three different types of reference evaporation from data of KNMI meteo stations:\n",
    "\n",
    "- Penman\n",
    "- Hargreaves\n",
    "- Makkink (in the same way as KNMI)\n",
    "\n",
    "These three types of reference evaporation are calculated the same way as described by [Allen et al. 1990](https://www.fao.org/3/x0490E/x0490e07.htm#solar%20radiation) and [STOWA rapport](https://edepot.wur.nl/163482). Be aware that the last report is written in Dutch and contains errors in the units.\n",
    "\n",
    "The following variables from the KNMI are used for each reference evaporation type:\n",
    "\n",
    "- Penman: average (TG), minimum (TN) and maximum (TX) temperature, de global radiation (Q), de windspeed (FG) en de relative humidity (PG).\n",
    "- Hargreaves: average (TG), minimum (TN) and maximum (TX) temperature\n",
    "- Makkink: average temperature (TG) and global radiation (Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison Makkink\n",
    "\n",
    "Lets compare Hypdropandas Makkink verdamping evaporation with the EV24 Makkink verdamping of the KNMI. When Hydropandas Makkink evaporation is rounded (on 4 decimals), the estimate is the same as for the KNMI estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev24 = hpd.EvaporationObs.from_knmi(\n",
    "    stn=260, start=\"2022-1-1\"\n",
    ")  # et_type='EV24' by default\n",
    "makk = hpd.EvaporationObs.from_knmi(meteo_var=\"makkink\", stn=260, start=\"2022-1-1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(2, figsize=(11, 4))\n",
    "ax[0].plot(ev24, label=ev24.name)\n",
    "ax[0].plot(makk.round(4), label=makk.name)\n",
    "ax[0].set_ylabel(\"E [m/d]\")\n",
    "ax[0].set_title(\"Makkink evaporation\")\n",
    "ax[0].legend()\n",
    "ax[1].plot(ev24[\"EV24\"] - makk[\"makkink\"].round(4))\n",
    "ax[1].set_title(\"Difference Makkink KNMI and Hydropandas\")\n",
    "f.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison Penman, Makkink en Hargreaves\n",
    "On average Penman gives a higher estimate for reference evaporation than Makkink (~0.55mm). This can be explained by the fact that Penman takes into account windspeed and Makkink ignores this proces. Hargreaves is a very simple way of estimation the evaporation, only taking into account temperature and extraterrestial radiation. Therefore it gives a lower estimate for the reference evporatoin compared to the two other methods (~-0.35mm wrt Makkink)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "penm = hpd.EvaporationObs.from_knmi(\n",
    "    meteo_var=\"penman\", stn=260, start=\"2022-1-1\"\n",
    ").squeeze()\n",
    "harg = hpd.EvaporationObs.from_knmi(\n",
    "    meteo_var=\"hargreaves\", stn=260, start=\"2022-1-1\"\n",
    ").squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(11, 4))\n",
    "ax.plot(makk, label=makk.name)\n",
    "ax.plot(penm, label=penm.name)\n",
    "ax.plot(harg, label=harg.name)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Spatial interpolate (Makkink) Evaporation?](#top)<a id=EvapInterp></a>\n",
    "\n",
    "Does interpolation matter? There are ways to interpolate evaporation datasets. However currently the nearest station is always used in HydroPandas' methods. Does this give different results? First lets look spatially."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get all stations where EV24 is measured"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stns = knmi.get_stations(meteo_var=\"EV24\").sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collect all EV24 data ever measured by KNMI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = \"1950-01-01\"\n",
    "tmax = \"2022-04-11\"\n",
    "\n",
    "# empty dataframe\n",
    "df = pd.DataFrame(\n",
    "    columns=stns.index\n",
    ")  # index=pd.date_range(start=tmin, end=tmax, freq='H')\n",
    "\n",
    "for stn in tqdm(stns.index):\n",
    "    df_stn = hpd.MeteoObs.from_knmi(\n",
    "        meteo_var=\"EV24\", stn=stn, fill_missing_obs=False, start=tmin, end=tmax\n",
    "    )\n",
    "    df[stn] = df_stn\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the KNMI, thin plate spline is the best way to interpolate Makkink evaporation. Thats also how they provide the gridded Makkink evaporation : \n",
    "\n",
    "- [Evaporation Dataset - gridded daily Makkink evaporation for the Netherlands](https://dataplatform.knmi.nl/dataset/ev24-2)\n",
    "- [Interpolation of Makkink evaporation in the\n",
    "Netherlands - Paul Hiemstra and Raymond Sluiter (2011)](https://cdn.knmi.nl/knmi/pdf/bibliotheek/knmipubTR/TR327.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy = stns.loc[df.columns, [\"x\", \"y\"]]\n",
    "\n",
    "for idx in tqdm(df.index[0 : len(df) : 2000]):\n",
    "    # get all stations with values for this date\n",
    "    val = df.loc[idx].dropna() * 1000  # mm\n",
    "    # get stations for this date\n",
    "    coor = xy.loc[val.index].to_numpy()\n",
    "    if (\n",
    "        len(val) < 3\n",
    "    ):  # if there are less than 3 stations, thin plate spline does not work\n",
    "        # options: linear, multiquadric, gaussian,\n",
    "        kernel = \"linear\"\n",
    "\n",
    "    else:\n",
    "        kernel = \"thin_plate_spline\"\n",
    "        # options:\n",
    "        # 'inverse_quadratic', 'linear', 'multiquadric', 'gaussian',\n",
    "        # 'inverse_multiquadric', 'cubic', 'quintic', 'thin_plate_spline'\n",
    "\n",
    "    # create an scipy interpolator\n",
    "    rbf = RBFInterpolator(coor, val.to_numpy(), epsilon=1, kernel=kernel)\n",
    "\n",
    "    nea = NearestNDInterpolator(coor, val.to_numpy())\n",
    "\n",
    "    # interpolate on grid of the Netherlands\n",
    "    grid = np.mgrid[10000:280000:100j, 300000:620000:100j]\n",
    "    grid2 = grid.reshape(2, -1).T  # interpolator only takes array [[x0, y0],\n",
    "    #  [x1, y1]]\n",
    "    val_rbf = rbf.__call__(grid2).reshape(100, 100)\n",
    "    val_nea = nea.__call__(grid2).reshape(100, 100)\n",
    "\n",
    "    # create figure\n",
    "    fig, ax = plt.subplot_mosaic(\"AAAABBBBC\", figsize=(10, 5.925))\n",
    "    fig.suptitle(f\"Makkink Evaporation [mm] {idx.date()}\", y=0.95)\n",
    "    vmin = 0\n",
    "    vmax = 5\n",
    "\n",
    "    ax[\"A\"].set_title(f\"Interpolation: {kernel}\")\n",
    "    ax[\"A\"].pcolormesh(grid[0], grid[1], val_rbf, vmin=vmin, vmax=vmax)\n",
    "    ax[\"B\"].set_title(\"Interpolation: Nearest\")\n",
    "    ax[\"B\"].pcolormesh(grid[0], grid[1], val_nea, vmin=vmin, vmax=vmax)\n",
    "    ax[\"A\"].scatter(*coor.T, c=val, s=100, ec=\"k\", vmin=vmin, vmax=vmax)\n",
    "    p = ax[\"B\"].scatter(*coor.T, c=val, s=100, ec=\"k\", vmin=vmin, vmax=vmax)\n",
    "    cb = fig.colorbar(p, cax=ax[\"C\"])\n",
    "    cb.set_label(\"[mm]\")\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same method is implemented in Hydropandas for an ObsCollection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sd = \"2022-01-01\"\n",
    "ed = \"2022-12-31\"\n",
    "oc = hpd.read_knmi(\n",
    "    stns=stns.index, starts=sd, ends=ed, meteo_vars=[\"EV24\"], raise_exceptions=False\n",
    ")\n",
    "oc_et = oc.interpolate(xy=[(100000, 330000)])\n",
    "eti = oc_et.iloc[0].obs\n",
    "eti"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "etn = hpd.MeteoObs.from_knmi(\n",
    "    xy=(100000, 330000), start=sd, end=ed, meteo_var=\"EV24\", fill_missing_obs=True\n",
    ")\n",
    "etn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, for one location the interpolation method is significantly slower. Lets see how the values compare for a time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True)\n",
    "eti.plot(ax=ax[0])\n",
    "etn.plot(ax=ax[0], linestyle=\"--\", color=\"C1\")\n",
    "ax[0].set_title(\"Comparison Interpolated and Nearest Makkink Evaporation\")\n",
    "ax[0].set_ylabel(\"[mm]\")\n",
    "ax[0].grid()\n",
    "ax[0].legend([\"Interpolated\", \"Nearest\"])\n",
    "\n",
    "(eti.squeeze() - etn[\"EV24\"].squeeze()).plot(ax=ax[1], color=\"C2\")\n",
    "ax[1].set_ylabel(\"[mm]\")\n",
    "ax[1].grid()\n",
    "ax[1].legend([\"Interpolated - Nearest\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolated evaporation can also be collected for multiple points (using x and y in a list of DataFrame) in an ObsCollection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oc_et = oc.interpolate(xy=[(100000, 320000), (100000, 330000), (100000, 340000)])\n",
    "oc_et"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpolation method is slow at first, but if collected for many different locations the time penalty is not that significant anymore. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}