{ "cells": [ { "cell_type": "markdown", "id": "87fd561f", "metadata": {}, "source": [ "# KNMI Drought Indices\n", "\n", "J.P.M. Witte, G.A.P.H. van den Eertwegh and P.J.J.F. Torfs (2025) - [Absolute Meteorological Drought Indices Validated Against Irrigation Amounts](https://doi.org/10.3390/w17071056)." ] }, { "cell_type": "markdown", "id": "5f4ec1ab", "metadata": {}, "source": [ "## Load packages" ] }, { "cell_type": "code", "execution_count": null, "id": "1a439ca7", "metadata": {}, "outputs": [], "source": [ "import matplotlib as mpl\n", "import pandas as pd\n", "\n", "from spei import knmi\n", "from spei.plot import deficit_knmi\n", "from spei.utils import group_yearly_df" ] }, { "cell_type": "markdown", "id": "ae00001a", "metadata": {}, "source": [ "## Get data\n", "Data from De Bilt (260) or P13 stations from 1960 till today" ] }, { "cell_type": "markdown", "id": "fba32ea8", "metadata": {}, "source": [ "### Most recent De Bilt data" ] }, { "cell_type": "code", "execution_count": null, "id": "131815ed", "metadata": {}, "outputs": [], "source": [ "# import hydropandas as hpd\n", "# prec = hpd.PrecipitationObs.from_knmi(\n", "# meteo_var=\"RH\",\n", "# stn=260,\n", "# startdate=pd.Timestamp(\"1960-01-01\"),\n", "# enddate=pd.Timestamp.today(),\n", "# )[\"RH\"].multiply(1e3)\n", "# prec.index = prec.index.normalize()\n", "# evap = hpd.EvaporationObs.from_knmi(\n", "# meteo_var=\"EV24\",\n", "# stn=260,\n", "# startdate=pd.Timestamp(\"1960-01-01\"),\n", "# enddate=pd.Timestamp.today(),\n", "# )[\"EV24\"].multiply(1e3)\n", "# evap.index = evap.index.normalize()\n", "# temp = hpd.MeteoObs.from_knmi(\n", "# meteo_var=\"TG\",\n", "# stn=260,\n", "# startdate=pd.Timestamp(\"1960-01-01\"),\n", "# enddate=pd.Timestamp.today(),\n", "# )[\"TG\"]\n", "# temp.index = temp.index.normalize()" ] }, { "cell_type": "markdown", "id": "72136d23", "metadata": {}, "source": [ "### KNMI stations data\n", "https://www.knmi.nl/kennis-en-datacentrum/achtergrond/achtergrondinformatie-klimaatdashboard\n", "\n", "De neerslagtekort klimaatdashboardgrafiek is alleen voor het landelijk gemiddelde beschikbaar, gebaseerd op:\n", "\n", "Voor 1906 t/m 2000: officiële reeks voor Nederland voor neerslagtekort: Dagelijks neerslagtekort NL (1 apr t/m 30 sep) op basis van Makkink verdamping De Bilt geschat uit zonneschijnduur minus 13 neerslagstations (P13) (c) KNMI, mei 2020, Jules Beersma: Climate Explorer \n", "\n", "Vanaf 2001: verdamping (gemiddelde van 13 automatische weerstations nabij 13 neerslagstations) minus de gemiddelde neerslag van 13 neerslagstations:\n", "De P13: het gemiddelde van de hoeveelheid neerslag op de volgende 13 KNMI-neerslagstations: De Bilt (550_N), De Kooy (25_N), Groningen (139_N), Heerde (328_N), Hoofddorp (438_N), Hoorn (222_N), Kerkwerve (737_N), Oudenbosch (828_N), Roermond (961_N), Ter Apel (144_N), West-Terschelling (11_N), Westdorpe (770_N) en Winterswijk (666_N).\n", "\n", "Het gemiddelde van de hoeveelheid verdamping (EV24) op 13 automatische weerstations van het KNMI nabij de 13 neerslagstations: De Bilt (260_H), De Kooy (235_H), Eelde (280_H), Heino (278_H), Schiphol (240_H), Berkhout (249_H), Vlissingen (310_H), Eindhoven (370_H), Ell (377_H), Nieuw Beerta (286_H), Hoorn Terschelling (251_H), Westdorpe (319_H) en Hupsel (283_H).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7401ab30", "metadata": {}, "outputs": [], "source": [ "# import hydropandas as hpd\n", "\n", "# P13 stations\n", "# p_stns = [\n", "# 550, # De Bilt\n", "# 25, # De Kooy\n", "# 139, # Groningen\n", "# 328, # Heerde\n", "# 438, # Hoofddorp\n", "# 222, # Hoorn\n", "# 737, # Kerkwerve\n", "# 828, # Oudenbosch\n", "# 961, # Roermond\n", "# 144, # Ter Apel\n", "# 11, # West-Terschelling\n", "# 770, # Westdorpe\n", "# 666, # Winterswijk\n", "# ]\n", "\n", "# # EV24-13 stations\n", "# ev_stns = [\n", "# 260, # De Bilt\n", "# 235, # De Kooy\n", "# 280, # Eelde\n", "# 278, # Heino\n", "# 240, # Schiphol\n", "# 249, # Berkhout\n", "# 310, # Vlissingen\n", "# 370, # Eindhoven\n", "# 377, # Ell\n", "# 286, # Nieuw Beerta\n", "# 251, # Hoorn Terschelling\n", "# 319, # Westdorpe\n", "# 283, # Hupsel\n", "# ]\n", "# oc_p = hpd.ObsCollection.from_knmi(\n", "# stns=p_stns,\n", "# starts=pd.Timestamp(\"1960-01-01\"),\n", "# ends=pd.Timestamp.today(),\n", "# meteo_vars=[\"RD\"],\n", "# )\n", "# oc_ev = hpd.ObsCollection.from_knmi(\n", "# stns=ev_stns,\n", "# starts=pd.Timestamp(\"1960-01-01\"),\n", "# ends=pd.Timestamp.today(),\n", "# meteo_vars=[\"EV24\", \"TG\"],\n", "# )\n", "\n", "# prec = pd.concat([o[\"RD\"] for o in oc_p[\"obs\"]], axis=1).mean(axis=1).multiply(1e3).rename(\"prec\")\n", "# prec.index = prec.index.normalize()\n", "# ev_data = pd.DataFrame({mv: pd.concat([o[mv] for o in gr[\"obs\"]], axis=1).mean(axis=1) for mv, gr in oc_ev.groupby(\"meteo_var\")})\n", "# ev_data.index = ev_data.index.normalize()\n", "# ev_data = ev_data.loc[prec.index] # align indices because prec stations less frequently reported\n", "# evap = ev_data[\"EV24\"].multiply(1e3).rename(\"evap\")\n", "# temp = ev_data[\"TG\"].rename(\"temp\")" ] }, { "cell_type": "markdown", "id": "24a7519c", "metadata": {}, "source": [ "### From a CSV file" ] }, { "cell_type": "code", "execution_count": null, "id": "528daaa6", "metadata": {}, "outputs": [], "source": [ "df = pd.read_csv(\"data/DEBILT.csv\", index_col=0, parse_dates=True)\n", "prec = df[\"Prec [m/d] 260_DEBILT\"].multiply(1e3).rename(\"prec\")\n", "evap = df[\"Evap [m/d] 260_DEBILT\"].multiply(1e3).rename(\"evap\")\n", "temp = df[\"Temp [C] 260_DEBILT\"].rename(\"temp\")" ] }, { "cell_type": "markdown", "id": "b268e457", "metadata": {}, "source": [ "## Calculate precipitation deficit" ] }, { "cell_type": "code", "execution_count": null, "id": "433a3701", "metadata": {}, "outputs": [], "source": [ "deficit = evap - prec\n", "## deficit period\n", "startdate = pd.Timestamp(\"2000-04-01\")\n", "enddate = pd.Timestamp(\"2000-09-30\")\n", "# calculate cumulative deficit\n", "cumdf = knmi.get_cumulative_deficit(\n", " deficit=deficit,\n", " startdate=startdate,\n", " enddate=enddate,\n", " allow_below_zero=False,\n", ")\n", "# plot deficit\n", "ax = cumdf.plot(figsize=(7.0, 6.0), cmap=\"cividis\")\n", "ax.legend(ncol=5, loc=(0, 1))\n", "ax.xaxis.set_major_locator(mpl.dates.MonthLocator())\n", "ax.xaxis.set_major_formatter(mpl.dates.DateFormatter(\"%B\"))\n", "ax.xaxis.set_ticks([], minor=True)\n", "ax.set_ylim(0.0)" ] }, { "cell_type": "markdown", "id": "b9a26454", "metadata": {}, "source": [ "## Precipitation deficit indices" ] }, { "cell_type": "code", "execution_count": null, "id": "a0a18964", "metadata": {}, "outputs": [], "source": [ "doct1 = knmi.deficit_oct1(deficit)\n", "doct1.to_frame().transpose()" ] }, { "cell_type": "code", "execution_count": null, "id": "ca3b85aa", "metadata": {}, "outputs": [], "source": [ "dmax = knmi.deficit_max(deficit)\n", "dmax.to_frame().transpose()" ] }, { "cell_type": "code", "execution_count": null, "id": "88fa3017", "metadata": {}, "outputs": [], "source": [ "diapr1 = knmi.deficit_apr1(deficit)\n", "diapr1.to_frame().transpose()" ] }, { "cell_type": "code", "execution_count": null, "id": "7fa91d05", "metadata": {}, "outputs": [], "source": [ "digdd = knmi.deficit_gdd(deficit, temp, threshold=440.0)\n", "digdd.to_frame().transpose()" ] }, { "cell_type": "code", "execution_count": null, "id": "66972df3", "metadata": {}, "outputs": [], "source": [ "diwet = knmi.deficit_wet(deficit)\n", "diwet.to_frame().transpose()" ] }, { "cell_type": "markdown", "id": "07085e76", "metadata": {}, "source": [ "## Compare to original KNMI data\n", "\n", "File obtained from https://climexp.knmi.nl/getindices.cgi?NPERYEAR=366&STATION=precipitationdeficit&TYPE=i&WMO=KNMIData/nt_nl&id=someone@somewhere" ] }, { "cell_type": "code", "execution_count": null, "id": "38d6f9af", "metadata": {}, "outputs": [], "source": [ "knmi_cumdf = group_yearly_df(\n", " pd.read_csv(\n", " \"data/neerslagtekort.txt\",\n", " skiprows=11,\n", " sep=\"\\t\",\n", " header=None,\n", " index_col=0,\n", " parse_dates=True,\n", " date_format=\"%Y%m%d\",\n", " )\n", " .dropna(how=\"all\", axis=1)\n", " .squeeze()\n", " .rename(\"KNMI\")\n", ")\n", "knmi_cumdf.index.name = \"\"\n", "ax = knmi_cumdf.plot(figsize=(7.0, 6.0), cmap=\"viridis\")\n", "ax.legend(ncol=5, loc=(0, 1))\n", "ax.xaxis.set_major_locator(mpl.dates.MonthLocator())\n", "ax.xaxis.set_major_formatter(mpl.dates.DateFormatter(\"%B\"))\n", "ax.xaxis.set_ticks([], minor=True)\n", "ax.set_ylim(0.0)" ] }, { "cell_type": "markdown", "id": "a5d7d7da", "metadata": {}, "source": [ "### KNMI plot\n", "\n", "From KNMI website the drought deficit is plotted as below:" ] }, { "cell_type": "markdown", "id": "a798d38d", "metadata": {}, "source": [ "![neerslagtekort](https://cdn.knmi.nl/knmi/map/page/klimatologie/grafieken/neerslagtekort/neerslagtekort.png)" ] }, { "cell_type": "markdown", "id": "22245717", "metadata": {}, "source": [ "#### With KNMI deficit data\n", "This plot can be reproduced (almost perfectly) as seen from the figure below.\n", "\n", "The calculation uses the average precipitation from 13 reference stations in the Netherlands (the so-called P13/EV24-13 stations) and the reference evaporation calculated based on sunshine duration in De Bilt (until 2001) or the global radiation near the P13 stations (from 2001 onwards). For the median and 5% driest years a rolling window is aplied. However, the size of this window is not documented anywhere." ] }, { "cell_type": "code", "execution_count": null, "id": "b06fa427", "metadata": {}, "outputs": [], "source": [ "ax = deficit_knmi(knmi_cumdf, window=28)\n", "ax.set_title(\"KNMI computed preciptiation deficit\")" ] }, { "cell_type": "markdown", "id": "e804afea", "metadata": {}, "source": [ "#### With own computed deficit (with downloaded knmi data)" ] }, { "cell_type": "code", "execution_count": null, "id": "eb4ea498", "metadata": {}, "outputs": [], "source": [ "ax = deficit_knmi(cumdf, window=0)\n", "ax.set_title(\"Downloaded measurements\")" ] } ], "metadata": { "kernelspec": { "display_name": "SPEI", "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.13.1" } }, "nbformat": 4, "nbformat_minor": 5 }