{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Drought Prediction with Time Series Modeling\n", "\n", "*Martin Vonk - 2022*\n", "\n", "This notebooks shows a quick calculation of the SPI, SPEI and SGI for De Bilt, in the Netherlands. The SGI is calculated using a [Pastas](https://github.com/pastas/pastas) time series model since the original time series is too short. The application of time series models for extrapolating groundwater time series is discussed in [Brakkee et al (2022)](https://hess.copernicus.org/articles/26/551/2022/hess-26-551-2022.html).\n", "\n", "## Required packages" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import pastas as ps\n", "import scipy.stats as scs\n", "\n", "import spei as si # si for standardized index\n", "\n", "print(si.show_versions())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import time series\n", "\n", "Time series are imported using the package hydropandas. Enddate is by default yesterday. The head time series is obtained from a Pastas test dataset." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import hydropandas as hpd\n", "\n", "# today = datetime.date.today()\n", "# yesterday = (today - datetime.timedelta(days=1)).strftime(\"%Y-%m-%d\")\n", "# prec = (\n", "# hpd.PrecipitationObs.from_knmi(\n", "# meteo_var=\"RH\", stn=260, startdate=\"1959-07-01\", enddate=yesterday\n", "# )\n", "# .multiply(1e3)\n", "# .squeeze()\n", "# )\n", "# prec.index = prec.index.normalize()\n", "# evap = (\n", "# hpd.EvaporationObs.from_knmi(\n", "# meteo_var=\"EV24\", stn=260, startdate=\"1959-07-01\", enddate=yesterday\n", "# )\n", "# .multiply(1e3)\n", "# .squeeze()\n", "# )\n", "# evap.index = evap.index.normalize()\n", "\n", "\n", "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", "head = df[\"Head [m] B32C0572_DEBILT\"].rename(\"B32C0572\").dropna()\n", "today = df.index[-1]\n", "yesterday = df.index[-2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate SPI and SPEI" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Accumulate time series on monthly basis\n", "spi1 = si.spi(prec, timescale=30, dist=scs.gamma, fit_freq=\"MS\")\n", "spei1 = si.spei((prec - evap), timescale=30, dist=scs.fisk, fit_freq=\"MS\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xlim = pd.to_datetime([\"2018-01-01\", df.index[-1]])\n", "\n", "fig, axs = plt.subplots(2, 1, figsize=(7.0, 5.5), sharex=True)\n", "si.plot.si(spi1, ybound=3.1, ax=axs[0], background=False, cmap=\"roma\")\n", "si.plot.si(spei1, ybound=3.1, ax=axs[1], background=False, cmap=\"roma\")\n", "[(x.grid(), x.set_xlim(xlim), x.set_ylabel(\"Z-Score\")) for x in axs]\n", "axs[0].set_title(\"Standardized Precipitation Index\")\n", "axs[1].set_title(\"Standardized Precipitation Evaporation Index\")\n", "fig.suptitle(\"Meteoroligical Drought-Indices De Bilt\")\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create time series model and simulate head " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ml = ps.Model(head)\n", "rm = ps.RechargeModel(\n", " prec, evap, ps.Exponential(), recharge=ps.rch.FlexModel(gw_uptake=True)\n", ")\n", "ml.add_stressmodel(rm)\n", "ml.solve(tmin=\"1970-07-01\", report=True)\n", "_ = ml.plots.results(figsize=(10.0, 8.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate SGI based on time series model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gws = ml.simulate(tmin=\"1990-07-01\", tmax=yesterday)\n", "sgi = si.sgi(gws, fit_freq=\"MS\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare three drought-indices (SPI, SPEI, SGI) in plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplot_mosaic(\n", " [[\"SPI\"], [\"SPEI\"], [\"SGI\"]], figsize=(6.5, 8), sharex=True\n", ")\n", "si.plot.si(spi1, ybound=3.5, ax=axs[\"SPI\"], add_category=False)\n", "si.plot.si(spei1, ybound=3.5, ax=axs[\"SPEI\"], add_category=False)\n", "si.plot.si(sgi, ybound=3.5, ax=axs[\"SGI\"], add_category=False)\n", "[(axs[x].grid(), axs[x].set(xlim=xlim, ylabel=\"Z-Score\")) for x in axs]\n", "axs[\"SPI\"].set_title(\"Standardized Precipitation Index 1\")\n", "axs[\"SPEI\"].set_title(\"Standardized Precipitation Evaporation Index 1\")\n", "axs[\"SGI\"].set_title(\"Standardized Groundwater Index\")\n", "fig.suptitle(\"Drought-Indices for De Bilt\", fontsize=14)\n", "fig.tight_layout()\n", "# fig.savefig('Drought_Index_Bilt.png', dpi=600, bbox_inches='tight')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare SPEI Kernel Density Estimate for one month" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax = si.plot.monthly_density(\n", " spi1, years=[today.year - 1, today.year], months=[today.month - 1]\n", ")\n", "ax.set_xlabel(\"Z-Score\")\n", "ax.set_title(\"SPEI\");" ] } ], "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.12.3" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }