07_resampling_and_aggregations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Grouping, resampling and aggregations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The hourly data over 23 years contain about ~800,000 observations. It is not viable to extract any knowledge by just looking at them ony by one. Thus an analysis of this sort of data almost always consists of proper **grouping** (selecting which rows we want to condense) and **aggregation** (applying an operation to produce a single or a handful of values from them)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Necessary import evil\n",
    "\n",
    "import pandas as pd\n",
    "from pandas import IndexSlice as idx, Grouper\n",
    "import pandera as pa\n",
    "from pandera.typing import Int16, DataFrame, Series, Category\n",
    "\n",
    "from weatherlyser.loader import load_chmi_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the saved data\n",
    "data = pd.read_parquet(\"./data/open_meteo_2000-2022.pq\")\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will skip this now as it takes a lot of time!\n",
    "# HistoricalWeatherDataFrame.validate(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Most meteorological reporting is done in local time\n",
    "data = data.tz_convert(level=\"time\", tz=\"Europe/Prague\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, we will be working with the \"best model\" for now. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model_data = data.loc[idx[:, \"best_match\"], :].reset_index(\"model\", drop=True)\n",
    "best_model_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weather description will become handy for our first aggregations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from weatherlyser.pa_models import WEATHER_CODES\n",
    "\n",
    "\n",
    "class WeatherCodeModel(pa.DataFrameModel):\n",
    "    weathercode: Int16 = pa.Field(ge=0, le=99, nullable=True)\n",
    "\n",
    "@pa.check_types\n",
    "# *** Annotations as a task\n",
    "def get_weather_description(df: DataFrame[WeatherCodeModel]) -> Series[Category]:\n",
    "    \"\"\"Converts the weather code column in a Pandas DataFrame\n",
    "    to a categorical column with the weather descriptions as the category levels.\n",
    "\n",
    "    Args:\n",
    "        df: A Pandas DataFrame with a column called \"weathercode\".\n",
    "\n",
    "    Returns:\n",
    "        A categorical series with a string representation of the code.\n",
    "    \"\"\"\n",
    "    # *** The mapping as a task\n",
    "    return df[\"weathercode\"].map(WEATHER_CODES).astype(\"category\")\n",
    "\n",
    "get_weather_description(best_model_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model_data = best_model_data.assign(\n",
    "    description=get_weather_description(best_model_data)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grouping and simple aggregations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped = best_model_data.groupby(\"description\")\n",
    "grouped"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is this object? Is this some numerical result yet? No, it is actually just the segmentation of the dataframe into **groups** of **rows** that will be handled by subsequent operation.\n",
    "\n",
    "The groupby mechanism ensures that the operations done on the sub-dataframes are recombined together with the grouping key(s) as index (levels of index):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a column and do a mean over it, within the categories\n",
    "grouped[\"rain\"].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiple aggregations over the same column\n",
    "grouped[\"rain\"].agg([\"min\", \"mean\", \"max\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With time series, we usually do not want to group over (all) values of column, but over some interval (typically of the index or index level). The mechanism is pretty much the same but the method is called `resample` and we need an extra argument, frequency.\n",
    "\n",
    "So in order to get each year's total precipitation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model_data.resample(\"1Y\")  # 1Y = every one year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This, in parallel to the groupby result is an object just holding the segmentations and any aggregation must be done on it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model_data.resample(\"1Y\")[\"precipitation\"].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will not be plotting much in this workshop but the simple `plot` method will still come handy to see what we calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model_data.resample(\"1Y\")[\"precipitation\"].sum().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Write a function `get_daily_temperature_stats` that finds the minimum, mean and maximum temperatures for each day. (Bonus: You can annotate it with a pandera model)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_daily_temperature_stats(df):\n",
    "   ...\n",
    "\n",
    "get_daily_temperature_stats(best_model_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, we want to combine grouping and resampling in one operation (although a nested grouping will probably work too). The **`Grouper`** class, used instead of strings as the `groupby` argument, allows to combine columns and normal/time index levels: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.groupby([Grouper(level=\"time\", freq=\"1D\"), Grouper(level=\"model\")])[\"precipitation\"].sum() #.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Have names for the aggregations: name=>(field, agg.method)\n",
    "data.groupby([Grouper(level=\"time\", freq=\"1D\"), Grouper(level=\"model\")]).agg(\n",
    "    min_pressure=(\"surface_pressure\", \"min\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Can you modify your get_daily_temperature_stats so that it works for the entire data, including model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_daily_temperature_stats2(df):\n",
    "\n",
    "\n",
    "\n",
    "    min_temp = (\"temp\", \"min\")\n",
    "    ...\n",
    "\n",
    "get_daily_temperature_stats2(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature engineering\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When looking at periodically repeating phenomena, we quite often do not want to \"just\" resample but look at aggregations over all the repeating intervals, such as when we are interested in (historical) mean temperature in May. The `resample` or `groupby` (even with Grouper) will not work straightaway.\n",
    "\n",
    "The easiest way is to construct helper columns with the time series feature (such as month number)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.index.get_level_values(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_temperature_per_month = best_model_data.assign(\n",
    "    month=lambda df: df.index.get_level_values(0).month\n",
    ").groupby(\"month\")[\"temperature_2m\"].mean()\n",
    "\n",
    "mean_temperature_per_month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_temperature_per_month.plot(kind=\"bar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise** For each calendar month (i.e. January, ...), find the typical (mean) hourly evolution of temperature over the course of the day. For that, you can use the attributes of the index and the method `unstack` to turn a level of a hierarchical index into columns (for each month)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def daily_monthly_pattern(df, column=\"temperature\"):\n",
    "    result = ...\n",
    "    return result.unstack(\"month\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "daily_monthly_pattern(best_model_data).plot()\n",
    "daily_monthly_pattern(best_model_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Complex exercise**: For each year, find how many minimum and maximum temperature days it has (the media tell us the weather is becoming more and more extreme!)\n",
    "\n",
    "We will do this in multiple steps:\n",
    "\n",
    "* Construct features \"year\", \"month\" and \"day\".\n",
    "* Group over those and find the extremes, use `idxmax`, `idxmin` methods.\n",
    "* Extract the year from both series (beware, that you need to use the `.dt` accessor for series datetime methods)\n",
    "* Use value_counts or groupby to count the year occurences\n",
    "\n",
    "Optional parts:\n",
    "* Create a DataFrameModel checking the output: int as index (should be consecutive years somewhere around 2000), two columns (min_count, max_count?) as positive integes\n",
    "* Create a hypothesis test showing that the total number of extremes does not exceed 366 over a dataframe with more than one-year worth of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_year_extremes(temperature_data: pd.DataFrame) -> pd.DataFrame:\n",
    "    # Create the year, month, day features\n",
    "    temperature_With_features = ...\n",
    "    \n",
    "    # For each month / day, get the extreme values\n",
    "    extreme_dates = ...\n",
    "\n",
    "    # Extract the years component\n",
    "    extreme_years = ...\n",
    "\n",
    "    # Do the final counting\n",
    "    extreme_counts = ...\n",
    "\n",
    "    return extreme_counts\n",
    "\n",
    "\n",
    "find_year_extremes(temperature_data=get_daily_temperature_stats(best_model_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# On CHMI data\n",
    "\n",
    "chmi_data = load_chmi_data()\n",
    "chmi_data_renamed = chmi_data.rename(columns={\"minimum_temperature\": \"min\", \"maximum_temperature\": \"max\"})\n",
    "find_year_extremes(chmi_data_renamed)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}