{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "\n",
    "# Heatwave detection\n",
    "\n",
    "According to the Intergovernmental Panel on Climate Change (IPCC), a heatwave can be defined as \"a period of abnormally hot weather, often defined with reference to a relative temperature threshold, lasting from two days to months\" {cite}`ipcc2022annex`. Such a definition is rather vague and can be easily adapted to depending on the local context. For instance, MeteoSwiss distinguishes between four level of heat warnings {cite}`burgstall2021meteoswiss` depending on both the period duration and temperature thresholds.\n",
    "\n",
    "This notebook shows how meteora can be used to detect heatwave periods, providing the user with flexibility to adapt to local heatwave definitions. To that end, we will use the [Swiss Federal Office of Meteorology and Climatology (MeteoSwiss)](https://www.meteoswiss.admin.ch) dataset. Let us start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "\n",
    "import contextily as cx\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "from meteora import clients, utils\n",
    "\n",
    "figwidth, figheight = plt.rcParams[\"figure.figsize\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "We will get the MeteoSwiss data for the June-July-August (JJA) periods of 2022 to 2024 for an area surrounding the city of Zurich (Switzerland):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "region = \"Zurich, Switzerland\"\n",
    "# select study period\n",
    "start_year = 2022\n",
    "end_year = 2024\n",
    "# months to consider when querying the MeteoSwiss data\n",
    "start_month = 6\n",
    "end_month = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "We can first plot the station locations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "client = clients.MeteoSwissClient(region)\n",
    "ax = client.stations_gdf.plot()\n",
    "cx.add_basemap(ax, crs=client.stations_gdf.crs, attribution=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "*(C) OpenStreetMap contributors, Tiles style by Humanitarian OpenStreetMap Team hosted by OpenStreetMap France*\n",
    "\n",
    "Then, we will get the JJA time series of temperature measurements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_df = pd.concat(\n",
    "    [\n",
    "        client.get_ts_df(\n",
    "            \"temperature\",\n",
    "            dt.date(year, start_month, 1),\n",
    "            # get latest moment of the latest day of the month\n",
    "            # ACHTUNG: this won't work if `end_month` is 12 (see commented code below)\n",
    "            dt.datetime.combine(\n",
    "                dt.date(year, end_month + 1, 1) - dt.timedelta(days=1), dt.time.max\n",
    "            ),\n",
    "            # dt.date(year + int(end_month / 12), (end_month % 12) + 1, 1)\n",
    "            # - dt.timedelta(days=1),\n",
    "        )\n",
    "        for year in range(start_year, end_year + 1)\n",
    "    ],\n",
    "    axis=\"rows\",\n",
    ")\n",
    "# convert to wide format\n",
    "ts_df = utils.long_to_wide(ts_df)\n",
    "ts_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "According to the MeteoSwiss heat warning system, a *level 4* warning is issued when the daily mean temperature is equal or greater than 27°C for at least 3 days. In order to detect a period with such characteristics, we can use the `meteora.utils.get_heatwave_periods` function, which takes a time series data frame as the only positional argument, followed by a series of keyword-only arguments. For this example, we need to set:\n",
    "- `heatwave_t_threshold` as the temperature threshold (the default is 25)\n",
    "- `heatwave_n_consecutive_days` as the number of consecutive days (the default is 3)\n",
    "- `station_agg_func` as the function to aggregate each station's daily time series (the default is \"mean\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_periods = utils.get_heatwave_periods(\n",
    "    ts_df,\n",
    "    heatwave_t_threshold=27,\n",
    "    heatwave_n_consecutive_days=3,\n",
    "    station_agg_func=\"mean\",\n",
    ")\n",
    "heatwave_periods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "As we can see, there is no heatwave fullfilling these characteristics for our study area during the specified period. Let us now try to find periods of *level 3* warnings, i.e., when the mean daily temperature is equal or greater than 25°C for at least 3 days:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_periods = utils.get_heatwave_periods(\n",
    "    ts_df,\n",
    "    heatwave_t_threshold=25,\n",
    "    heatwave_n_consecutive_days=3,\n",
    "    station_agg_func=\"mean\",\n",
    ")\n",
    "heatwave_periods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "We get a single heatwave that spans from the 19th to the 24th of August 2023. We can also use this information to extract the time series of this period:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_start, heatwave_end = heatwave_periods[0]\n",
    "ts_df.loc[heatwave_start:heatwave_end].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "Since there are three stations in our study area, we have three daily means to compare with `heatwave_t_threshold`. By default, the `meteora.utils.get_heatwave_periods` function will compare the threshold with the mean among all stations of `ts_df`, however we can change this behaviour through the `inter_station_agg_func` keyword argument. For instance, we can set to `\"max\"` so that a day is considered as a potential part of a heatwave when (at least) one of the stations is over the temperature threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_periods = utils.get_heatwave_periods(\n",
    "    ts_df,\n",
    "    heatwave_t_threshold=25,\n",
    "    heatwave_n_consecutive_days=3,\n",
    "    station_agg_func=\"mean\",\n",
    "    inter_station_agg_func=\"max\",\n",
    ")\n",
    "heatwave_periods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "When using such a criteria we obtain two heatwaves instead of one. The function intents to offer flexibility so that the user can easily adapt it to different heatwave definitions. For instance, we could also change the `station_agg_func` to `\"max\"` to compare the temperature threshold with the maximum daily temperature of each station, instead of the mean. However, the mean temperature is often most correlated with epidemiological impacts of heat stress because it also reflects night-time temperatures {cite}`burgstall2021meteoswiss`.\n",
    "\n",
    "Instead of getting the periods, we can use the `meteora.utils.get_heatwave_ts_df` function to get their corresponding time series data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_ts_df = utils.get_heatwave_ts_df(ts_df, heatwave_periods=heatwave_periods)\n",
    "heatwave_ts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "We now get a multi-indexed time series data frame where the outermost index is the heatwave identifier and the innermost index is the time. We can use this data frame, e.g., to plot the time series of each heatwave separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_gb = heatwave_ts_df.groupby(\"heatwave\")\n",
    "fig, axes = plt.subplots(\n",
    "    nrows=1, ncols=len(heatwave_gb), figsize=(len(heatwave_gb) * figwidth, figheight)\n",
    ")\n",
    "for (heatwave, _ts_df), ax in zip(heatwave_gb, axes):\n",
    "    _ts_df.droplevel(\"heatwave\").plot(ax=ax)\n",
    "    ax.set_title(f\"Heatwave {heatwave}\")\n",
    "    ax.set_ylabel(\"Temperature (°C)\")\n",
    "    ax.set_xlabel(\"Time\")\n",
    "    # rotate x-axis labels\n",
    "    ax.tick_params(axis=\"x\", rotation=45)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "We can also use the `meteora.utils.get_heatwave_ts_df` to get the time series data frame for the heatwave periods directly, without having to first get the periods. To that end, instead of providing the `heatwave_periods` argument, we can use the same keyword arguments as in the `meteora.utils.get_heatwave_periods` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatwave_ts_df = utils.get_heatwave_ts_df(\n",
    "    ts_df,\n",
    "    heatwave_t_threshold=25,\n",
    "    heatwave_n_consecutive_days=3,\n",
    "    station_agg_func=\"mean\",\n",
    "    inter_station_agg_func=\"max\",\n",
    ")\n",
    "heatwave_ts_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "```{bibliography}\n",
    ":filter: docname in docnames\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "tags,-all"
  },
  "kernelspec": {
   "display_name": "Python (Pixi)",
   "language": "python",
   "name": "pixi-kernel-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.11"
  },
  "pixi-kernel": {
   "environment": "doc"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}