{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G8aCTFO38ACI"
      },
      "source": [
        "# Introduction\n",
        "\n",
        "Italy compresses alpine, continental, and Mediterranean climates into one narrow peninsula. That makes it a good place to study how warming and rainfall shifts play out across very different local environments.\n",
        "\n",
        "In this post I use the updated ISTAT meteoclimatic workbook for **provincial capitals** to answer a few simple questions:\n",
        "\n",
        "- How quickly are Italian provincial capitals warming?\n",
        "- Which places are hottest on average?\n",
        "- Is precipitation moving in the same direction?\n",
        "- What did 2023 look like relative to recent and historical baselines?\n",
        "\n",
        "The source workbook is rich, but awkwardly structured. As often happens with official statistics, the hardest part is not plotting the charts — it is turning a human-readable Excel layout into tidy data.\n",
        "\n",
        "## A Note for Aspiring Data Scientists\n",
        "\n",
        "One of the most powerful skills a data scientist can develop is the ability to turn raw data into clear, actionable insights. National statistical bodies like **ISTAT** in Italy, **Eurostat**, or the **US Census Bureau** provide a wealth of free, high-quality datasets to experiment and gain insights from.\n",
        "\n",
        "As you'll see in this analysis, the technical barrier is actually quite low. Once you learn how to handle structured workbooks and perform basic Exploratory Data Analysis (EDA), you can uncover stories about the fields you are passionate about—from climate change to economics or social trends. I hope this notebook serves as a roadmap for your own discovery process.\n",
        "\n",
        "*This article rebuilds the original notebook-style analysis using the updated ISTAT workbook for provincial capitals, covering annual temperature and precipitation from 2006 to 2023, plus anomaly tables for 2023.*\n",
        "\n",
        "If you like this work, please give a star to the [related Github repo by following this link](https://www.github.com/feddernico/temperature-italia)\n",
        "\n",
        "# Loading the workbook\n",
        "\n",
        "The key sheets for the core analysis are:\n",
        "\n",
        "- `Tavola_1`: annual mean temperature by provincial capital, 2006–2023\n",
        "- `Tavola_2`: annual precipitation by provincial capital, 2006–2023\n",
        "- `Tavola_3`: 2023 temperature anomalies versus 2006–2015 and climatological normals"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1YysTKJkWgIh",
        "outputId": "6cf42d37-1902-49f3-f055-4850555cc1ab"
      },
      "outputs": [],
      "source": [
        "!wget https://www.istat.it/wp-content/uploads/2024/06/Tavole_Dati-Meteoclimatici_Capoluoghi-Provincia_Anno-2022-e-serie-2006-2022_DCAT.xlsx\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WaB6aPluXgpJ",
        "outputId": "e1daa87d-ed54-44d8-da99-caf67cfca0db"
      },
      "outputs": [],
      "source": [
        "import matplotlib\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import plotly.express as px\n",
        "import plotly.graph_objects as go\n",
        "from pathlib import Path\n",
        "\n",
        "FILE_PATH = (\n",
        "    \"Tavole-_Dati-Meteoclimatici_Capoluoghi-Provincia_Anno-2023-e-serie-2006-2023.xlsx\"\n",
        ")\n",
        "\n",
        "# --- Blog design system palette (see design-system.html) -----------------\n",
        "BLUE = \"#2780e3\"\n",
        "CHARCOAL = \"#343a40\"\n",
        "INK = \"#1e1e1e\"\n",
        "TEXT = \"#343a40\"\n",
        "MUTED = \"rgba(52, 58, 64, 0.72)\"\n",
        "MUTED_HEX = \"#6c757d\"\n",
        "LINE = \"#dee2e6\"\n",
        "SOFT = \"#f8f9fa\"\n",
        "SOFT_BLUE = \"#d4e6f9\"\n",
        "CODE = \"#7d12ba\"\n",
        "\n",
        "matplotlib.rcParams.update(\n",
        "    {\n",
        "        \"font.family\": \"sans-serif\",\n",
        "        \"font.size\": 12,\n",
        "        \"text.color\": TEXT,\n",
        "        \"axes.labelcolor\": TEXT,\n",
        "        \"axes.edgecolor\": LINE,\n",
        "        \"axes.titlecolor\": TEXT,\n",
        "        \"axes.titleweight\": \"normal\",\n",
        "        \"axes.titlesize\": 16,\n",
        "        \"axes.spines.top\": False,\n",
        "        \"axes.spines.right\": False,\n",
        "        \"axes.grid\": True,\n",
        "        \"grid.color\": LINE,\n",
        "        \"grid.linewidth\": 0.8,\n",
        "        \"grid.linestyle\": \"-\",\n",
        "        \"axes.facecolor\": \"#ffffff\",\n",
        "        \"figure.facecolor\": \"#ffffff\",\n",
        "        \"xtick.color\": MUTED_HEX,\n",
        "        \"ytick.color\": MUTED_HEX,\n",
        "        \"axes.linewidth\": 1.0,\n",
        "        \"legend.frameon\": False,\n",
        "        \"legend.fontsize\": 11,\n",
        "        \"figure.dpi\": 130,\n",
        "    }\n",
        ")\n",
        "\n",
        "xls = pd.ExcelFile(FILE_PATH)\n",
        "xls.sheet_names"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gmJJEDtCX4IU"
      },
      "source": [
        "# From wide to tidy\n",
        "\n",
        "The annual tables are stored in **wide** form: one row per city and one column per year. That format is convenient for visual inspection, but awkward for analysis. A small parser lets us melt each sheet into long form while handling the source file’s missing-value sentinel (`....`)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 298
        },
        "id": "Yhb2xvwgXwsx",
        "outputId": "0cf345ed-9462-4752-993b-8a7e6c1c3faa"
      },
      "outputs": [],
      "source": [
        "def parse_table(sheet_name, value_name):\n",
        "    df = pd.read_excel(FILE_PATH, sheet_name=sheet_name, header=None)\n",
        "\n",
        "    years = [int(x) for x in df.iloc[3, 1:19].tolist()]\n",
        "    data = df.iloc[4:, :19].copy()\n",
        "    data.columns = [\"Province\"] + years\n",
        "\n",
        "    data = data[data[\"Province\"].notna()]\n",
        "    data = data[\n",
        "        ~data[\"Province\"]\n",
        "        .astype(str)\n",
        "        .str.contains(\n",
        "            r\"Fonte|Totale|Nota|Nelle tavole|^\\s*$\",\n",
        "            regex=True,\n",
        "            na=False,\n",
        "        )\n",
        "    ]\n",
        "\n",
        "    long = data.melt(id_vars=\"Province\", var_name=\"Year\", value_name=value_name)\n",
        "    long[\"Province\"] = long[\"Province\"].astype(str).str.strip()\n",
        "    long[\"Year\"] = pd.to_numeric(long[\"Year\"], errors=\"coerce\")\n",
        "    long[value_name] = pd.to_numeric(\n",
        "        long[value_name].replace(\"....\", np.nan),\n",
        "        errors=\"coerce\",\n",
        "    )\n",
        "    return long\n",
        "\n",
        "\n",
        "temperature_data = parse_table(\"Tavola_1\", \"Temperature\")\n",
        "precipitation_data = parse_table(\"Tavola_2\", \"Precipitation\")\n",
        "\n",
        "italy_weather_data = temperature_data.merge(\n",
        "    precipitation_data,\n",
        "    on=[\"Province\", \"Year\"],\n",
        "    how=\"inner\",\n",
        ")\n",
        "\n",
        "italy_weather_data.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "f28ed85f"
      },
      "outputs": [],
      "source": [
        "import warnings\n",
        "warnings.simplefilter(action='ignore', category=FutureWarning)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ELRuxt_wX8sA"
      },
      "source": [
        "A quick check confirms that the cleaned dataset spans the full 2006–2023 range."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 143
        },
        "id": "YvTBQan6X6Fu",
        "outputId": "25266a11-35a2-45df-f587-9e9fa0b9957e"
      },
      "outputs": [],
      "source": [
        "italy_weather_data.agg(\n",
        "    provinces=(\"Province\", \"nunique\"),\n",
        "    start_year=(\"Year\", \"min\"),\n",
        "    end_year=(\"Year\", \"max\"),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6xDuA3M0X_uy"
      },
      "source": [
        "# National trends\n",
        "\n",
        "A natural first view is the yearly mean across all provincial capitals. This compresses local variation into a single national summary and makes the direction of change much easier to see."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "u88lhzdiX-es",
        "outputId": "aed31b7b-4cb6-4376-e540-e47a4db912ec"
      },
      "outputs": [],
      "source": [
        "def linear_slope(df, y_col):\n",
        "    clean = df.dropna(subset=[y_col]).copy()\n",
        "    x = clean[\"Year\"].to_numpy(dtype=float)\n",
        "    y = clean[y_col].to_numpy(dtype=float)\n",
        "    if len(x) < 2:\n",
        "        return np.nan\n",
        "    return float(np.polyfit(x, y, 1)[0])\n",
        "\n",
        "\n",
        "italy_temperature_trend = (\n",
        "    italy_weather_data.groupby(\"Year\", as_index=False)[\"Temperature\"].mean().dropna()\n",
        ")\n",
        "\n",
        "temp_coef = np.polyfit(\n",
        "    italy_temperature_trend[\"Year\"].to_numpy(dtype=float),\n",
        "    italy_temperature_trend[\"Temperature\"].to_numpy(dtype=float),\n",
        "    1,\n",
        ")\n",
        "italy_temperature_trend[\"Trend\"] = (\n",
        "    temp_coef[0] * italy_temperature_trend[\"Year\"] + temp_coef[1]\n",
        ")\n",
        "\n",
        "national_temp_slope = linear_slope(italy_temperature_trend, \"Temperature\")\n",
        "national_temp_change = national_temp_slope * (2023 - 2006)\n",
        "\n",
        "national_temp_slope, national_temp_change"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WBBypm9_YDTb"
      },
      "source": [
        "cross provincial capitals, the fitted warming rate is about **0.055 °C per year**, or roughly **0.9–1.0 °C over the full window**. That does not mean every city warms at exactly the same pace, but it sets the national backdrop for the rest of the analysis.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 398
        },
        "id": "wnUEr_NiYBbD",
        "outputId": "155492df-f6de-4fc7-ec17-19bf054c27e8"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.ticker import MaxNLocator\n",
        "\n",
        "def parse_table(sheet_name, value_name):\n",
        "    df = pd.read_excel(FILE_PATH, sheet_name=sheet_name, header=None)\n",
        "    years = [int(x) for x in df.iloc[3, 1:19].tolist()]\n",
        "    data = df.iloc[4:, :19].copy()\n",
        "    data.columns = [\"Province\"] + years\n",
        "    data = data[data[\"Province\"].notna()]\n",
        "    data[\"Province\"] = data[\"Province\"].astype(str).str.strip()\n",
        "    data = data[~data[\"Province\"].str.contains(r\"Fonte|Totale|Nota|Nelle tavole|^$\", regex=True, na=False)]\n",
        "    long = data.melt(id_vars=\"Province\", var_name=\"Year\", value_name=value_name)\n",
        "    long[\"Year\"] = pd.to_numeric(long[\"Year\"], errors=\"coerce\")\n",
        "    long[value_name] = pd.to_numeric(long[value_name].replace(\"....\", np.nan), errors=\"coerce\")\n",
        "    return long\n",
        "\n",
        "italy_mean = italy_weather_data.groupby(\"Year\", as_index=False)[\"Temperature\"].mean().dropna()\n",
        "coef = np.polyfit(italy_mean[\"Year\"], italy_mean[\"Temperature\"], 1)\n",
        "italy_mean[\"Trend\"] = coef[0] * italy_mean[\"Year\"] + coef[1]\n",
        "\n",
        "BG = \"#f7f6f2\"; PANEL = \"#f9f8f5\"; TEXT = \"#28251d\"; MUTED = \"#7a7974\"; GRID = \"#dcd9d5\"; ACCENT = \"#01696f\"\n",
        "\n",
        "plt.rcParams.update({\n",
        "    \"figure.facecolor\": BG, \"axes.facecolor\": PANEL, \"axes.edgecolor\": GRID, \"axes.labelcolor\": TEXT,\n",
        "    \"axes.titlecolor\": TEXT, \"xtick.color\": MUTED, \"ytick.color\": MUTED, \"text.color\": TEXT,\n",
        "    \"font.size\": 11, \"axes.titlesize\": 18, \"axes.titleweight\": \"bold\", \"legend.frameon\": False\n",
        "})\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(10.5, 5.8))\n",
        "ax.plot(italy_mean[\"Year\"], italy_mean[\"Temperature\"], color=ACCENT, linewidth=2.8, marker=\"o\", markersize=5.5, label=\"Annual mean\", zorder=3)\n",
        "ax.plot(italy_mean[\"Year\"], italy_mean[\"Trend\"], color=TEXT, linewidth=1.8, linestyle=(0, (4, 4)), alpha=0.7, label=\"Linear trend\", zorder=2)\n",
        "\n",
        "# Set title with padding and place subtitle below it\n",
        "ax.set_title(\"Average annual temperature across Italian provincial capitals\", loc=\"left\", pad=38)\n",
        "ax.text(0, 1.02, \"Mean annual temperature by provincial capital, averaged nationally, 2006–2023\", transform=ax.transAxes, fontsize=10.5, color=MUTED)\n",
        "\n",
        "ax.set_xlabel(\"Year\"); ax.set_ylabel(\"Temperature (°C)\")\n",
        "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n",
        "ax.set_xticks(italy_mean[\"Year\"][::2])\n",
        "ax.grid(axis=\"y\", color=GRID, linewidth=0.8, alpha=0.85)\n",
        "ax.grid(axis=\"x\", visible=False)\n",
        "for spine in [\"top\", \"right\"]: ax.spines[spine].set_visible(False)\n",
        "ax.spines[\"left\"].set_color(GRID); ax.spines[\"bottom\"].set_color(GRID)\n",
        "ax.legend(loc=\"upper left\", ncol=2)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6EFVfApfm6Lv"
      },
      "source": [
        "Precipitation is much noisier. The annual mean can still be summarized with a simple line, but the year-to-year swings are much larger than for temperature, so the fitted trend should be read more cautiously. The shaded band in the chart below is a **descriptive spread** — one standard deviation of the yearly national means — added only as a visual guide to how wildly precipitation jumps from year to year. It is *not* a confidence interval or any kind of uncertainty estimate around the trend."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Dehv9sNWm8W2",
        "outputId": "47d20840-c9e2-457d-8f4e-8b8db29a7f4b"
      },
      "outputs": [],
      "source": [
        "italy_precipitation_trend = (\n",
        "    italy_weather_data.groupby(\"Year\", as_index=False)[\"Precipitation\"]\n",
        "    .mean()\n",
        "    .dropna()\n",
        ")\n",
        "\n",
        "prec_coef = np.polyfit(\n",
        "    italy_precipitation_trend[\"Year\"].to_numpy(dtype=float),\n",
        "    italy_precipitation_trend[\"Precipitation\"].to_numpy(dtype=float),\n",
        "    1,\n",
        ")\n",
        "\n",
        "italy_precipitation_trend[\"Trend\"] = (\n",
        "    prec_coef[0] * italy_precipitation_trend[\"Year\"] + prec_coef[1]\n",
        ")\n",
        "\n",
        "national_prec_slope = linear_slope(\n",
        "    italy_precipitation_trend,\n",
        "    \"Precipitation\"\n",
        ")\n",
        "print(f\"The overall trend is {round(national_prec_slope, 2)}mm\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 398
        },
        "id": "6rkNRyyVntJP",
        "outputId": "39047554-56ac-4f5d-cd68-d5c370f24fbf"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(10.5, 5.8))\n",
        "\n",
        "# Descriptive spread: ±1 standard deviation of the yearly national means.\n",
        "# This is a visual guide to year-to-year variability, NOT a confidence\n",
        "# interval or an estimate of uncertainty around the trend line.\n",
        "prec_std = italy_precipitation_trend[\"Precipitation\"].std()\n",
        "italy_precipitation_trend[\"Lower\"] = italy_precipitation_trend[\"Precipitation\"] - prec_std\n",
        "italy_precipitation_trend[\"Upper\"] = italy_precipitation_trend[\"Precipitation\"] + prec_std\n",
        "\n",
        "# Shaded band\n",
        "ax.fill_between(italy_precipitation_trend[\"Year\"], italy_precipitation_trend[\"Lower\"], italy_precipitation_trend[\"Upper\"], facecolor=ACCENT, alpha=0.12, linewidth=0, zorder=1, label=\"±1σ (year-to-year spread)\")\n",
        "\n",
        "# Lines\n",
        "ax.plot(italy_precipitation_trend[\"Year\"], italy_precipitation_trend[\"Precipitation\"], color=ACCENT, linewidth=2.8, marker=\"o\", markersize=5.5, label=\"Annual mean\", zorder=3)\n",
        "ax.plot(italy_precipitation_trend[\"Year\"], italy_precipitation_trend[\"Trend\"], color=TEXT, linewidth=1.8, linestyle=(0, (4, 4)), alpha=0.7, label=\"Linear trend\", zorder=2)\n",
        "\n",
        "# Titles and Subtitles\n",
        "ax.set_title(\"Italy precipitation trend (2006–2023)\", loc=\"left\", pad=38)\n",
        "ax.text(0, 1.02, \"Mean annual precipitation across provincial capitals from ISTAT Tavola_2\", transform=ax.transAxes, fontsize=10.5, color=MUTED)\n",
        "\n",
        "ax.set_xlabel(\"Year\"); ax.set_ylabel(\"Rain (mm)\")\n",
        "ax.xaxis.set_major_locator(MaxNLocator(integer=True))\n",
        "ax.set_xticks(italy_precipitation_trend[\"Year\"][::2])\n",
        "ax.grid(axis=\"y\", color=GRID, linewidth=0.8, alpha=0.85)\n",
        "ax.grid(axis=\"x\", visible=False)\n",
        "for spine in [\"top\", \"right\"]: ax.spines[spine].set_visible(False)\n",
        "ax.spines[\"left\"].set_color(GRID); ax.spines[\"bottom\"].set_color(GRID)\n",
        "ax.legend(loc=\"lower right\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q4j29GxRYXLr"
      },
      "source": [
        "# The hottest places\n",
        "\n",
        "Averaging temperature across the full 2006–2023 period gives a simple \"warmest cities\" ranking. Unsurprisingly, the leaders cluster in Sicily and the far south.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "D5eny3mYYE-B",
        "outputId": "b216ff2c-cdca-4161-bb11-a03592fcc92a"
      },
      "outputs": [],
      "source": [
        "top_hottest = (\n",
        "    italy_weather_data.groupby(\"Province\", as_index=False)[\"Temperature\"]\n",
        "    .mean()\n",
        "    .sort_values(\"Temperature\", ascending=False)\n",
        "    .head(5)\n",
        ")\n",
        "\n",
        "top_hottest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t5BlxM4RYrAM"
      },
      "source": [
        "Being hot is not the same as warming quickly. To estimate warming rates, I fit a linear slope for each city after keeping only places with at least 15 valid annual observations, which reduces the influence of patchy reporting histories.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 419
        },
        "id": "xRLrLd7vYcLH",
        "outputId": "e1a7eda2-97f6-4eda-a7a2-e403a0cd46a5"
      },
      "outputs": [],
      "source": [
        "def slope_per_city(group, value_col):\n",
        "    clean = group.dropna(subset=[value_col])\n",
        "    x = clean[\"Year\"].to_numpy(dtype=float)\n",
        "    y = clean[value_col].to_numpy(dtype=float)\n",
        "    if len(x) < 2:\n",
        "        return np.nan\n",
        "    return round(float(np.polyfit(x, y, 1)[0]), 3)\n",
        "\n",
        "warming_base = italy_weather_data.groupby(\"Province\").filter(\n",
        "    lambda g: g[\"Temperature\"].notna().sum() >= 15\n",
        ")\n",
        "\n",
        "warming_rates = (\n",
        "    warming_base.groupby(\"Province\")\n",
        "    .apply(lambda g: slope_per_city(g, \"Temperature\"))\n",
        "    .rename(\"Trend (°C/year)\")\n",
        "    .reset_index()\n",
        "    .sort_values(\"Trend (°C/year)\", ascending=False)\n",
        ")\n",
        "warming_rates.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 446
        },
        "id": "85rvg0vuZdaR",
        "outputId": "841b5781-6b12-4255-fe5d-0aed8df0342d"
      },
      "outputs": [],
      "source": [
        "top_warming = warming_rates.head(10)\n",
        "fig, ax = plt.subplots(figsize=(10.5, 6.5))\n",
        "\n",
        "bars = ax.barh(top_warming[\"Province\"][::-1], top_warming[\"Trend (°C/year)\"][::-1], color=ACCENT, alpha=0.8, edgecolor=\"none\", height=0.65)\n",
        "\n",
        "# Subtitle below title\n",
        "ax.set_title(\"Fastest Warming Italian Provincial Capitals\", loc=\"left\", pad=38)\n",
        "ax.text(0, 1.02, \"Linear slope from annual mean temperature, cities with at least 15 valid years\", transform=ax.transAxes, fontsize=10.5, color=MUTED_HEX)\n",
        "\n",
        "ax.set_xlabel(\"Warming Trend (°C per year)\", color=MUTED_HEX, fontsize=11, labelpad=10)\n",
        "for spine in [\"top\", \"right\"]: ax.spines[spine].set_visible(False)\n",
        "ax.spines[\"left\"].set_color(GRID); ax.spines[\"bottom\"].set_color(GRID)\n",
        "ax.grid(axis=\"x\", linestyle=\"-\", color=GRID, alpha=0.8)\n",
        "ax.grid(axis=\"y\", visible=False)\n",
        "\n",
        "for bar in bars:\n",
        "    width = bar.get_width()\n",
        "    ax.text(width + 0.003, bar.get_y() + bar.get_height()/2, f\"+{width:.3f}\", va=\"center\", color=TEXT, fontsize=10, fontweight=\"bold\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lqg9ImcdedaW"
      },
      "source": [
        "# 2023 anomalies\n",
        "\n",
        "The long-run series are useful, but the workbook also includes a targeted anomaly table for 2023. `Tavola_3` compares 2023 against the 2006–2015 average and, where available, against the climatological normals for 1981–2010 and 1971–2000."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 262
        },
        "id": "G1lpoWY9efRn",
        "outputId": "4bb6cd16-0167-4ef3-fab2-7b7328f16e9d"
      },
      "outputs": [],
      "source": [
        "anomaly_sheet = pd.read_excel(FILE_PATH, sheet_name=\"Tavola_3\", header=None)\n",
        "\n",
        "anomaly_data = anomaly_sheet.iloc[4:114, :7].copy()\n",
        "anomaly_data.columns = [\n",
        "    \"Province\",\n",
        "    \"Diff_2023_vs_2006_2015\",\n",
        "    \"Mean_2006_2015\",\n",
        "    \"Anom_1981_2010\",\n",
        "    \"Clim_1981_2010\",\n",
        "    \"Anom_1971_2000\",\n",
        "    \"Clim_1971_2000\",\n",
        "]\n",
        "\n",
        "anomaly_data = anomaly_data[anomaly_data[\"Province\"].notna()]\n",
        "anomaly_data = anomaly_data[\n",
        "    ~anomaly_data[\"Province\"].astype(str).str.contains(\n",
        "        r\"Fonte|^\\s*$\",\n",
        "        regex=True,\n",
        "        na=False,\n",
        "    )\n",
        "]\n",
        "anomaly_data[\"Province\"] = anomaly_data[\"Province\"].astype(str).str.strip()\n",
        "\n",
        "for col in anomaly_data.columns[1:]:\n",
        "    anomaly_data[col] = pd.to_numeric(\n",
        "        anomaly_data[col].replace(\"....\", np.nan),\n",
        "        errors=\"coerce\",\n",
        "    )\n",
        "\n",
        "anomaly_data.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "31bb1f11"
      },
      "source": [
        "### Understanding the 1981–2010 Baseline\n",
        "\n",
        "In climate science, an **anomaly** is the difference between the observed temperature and a long-term average called a *climatological normal*. The World Meteorological Organization (WMO) uses 30-year periods to establish these normals.\n",
        "\n",
        "- **+0.0°C anomaly**: The year was exactly as warm as the 1981–2010 average.\n",
        "- **Positive anomaly (e.g., +2.66°C)**: The year was warmer than the historical baseline.\n",
        "\n",
        "When we see values exceeding +2.0°C in the table below, it indicates that those cities experienced a year significantly hotter than the conditions that characterized the end of the 20th century.\n",
        "\n",
        "A caveat on coverage: `Tavola_3` does not provide climatological-normal anomalies for every provincial capital. The sheet notes that those reference baselines are calculated only for **regional capitals** with long enough series, and some of the 2023-vs-2006–2015 differences are also missing where the recent series is incomplete. The ranking below is therefore **among places with available climatological comparisons, not all provincial capitals** — it shows which cities stand out within that comparable subset, not an exhaustive list of every warm place in Italy. Within that subset, 2023 stands out particularly strongly in cities such as **Perugia**, **Bologna**, **Milano**, and **Torino**."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "i4IAVXs6Y0_H",
        "outputId": "03efb0cb-4d49-4ae0-ab9f-77dc81faec74"
      },
      "outputs": [],
      "source": [
        "top_anomalies = (\n",
        "    anomaly_data.dropna(subset=[\"Anom_1981_2010\"])\n",
        "    .sort_values(\"Anom_1981_2010\", ascending=False)\n",
        "    [[\"Province\", \"Anom_1981_2010\"]]\n",
        "    .head(10)\n",
        ")\n",
        "\n",
        "top_anomalies"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 447
        },
        "id": "-DTBceQNfERr",
        "outputId": "6918477b-7024-4dc1-90ca-b8e65ddabcbf"
      },
      "outputs": [],
      "source": [
        "top_anomalies_plot = anomaly_data.dropna(subset=[\"Anom_1981_2010\"]).sort_values(\"Anom_1981_2010\", ascending=False).head(10)\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(10.5, 6.5))\n",
        "bars = ax.barh(top_anomalies_plot[\"Province\"][::-1], top_anomalies_plot[\"Anom_1981_2010\"][::-1], color=ACCENT, alpha=0.8, edgecolor=\"none\", height=0.65)\n",
        "\n",
        "# Subtitle below title\n",
        "ax.set_title(\"Largest 2023 temperature anomalies\", loc=\"left\", pad=38)\n",
        "ax.text(0, 1.02, \"Relative to 1981–2010 climatology from ISTAT Tavola_3\", transform=ax.transAxes, fontsize=10.5, color=MUTED_HEX)\n",
        "\n",
        "ax.set_xlabel(\"Anomaly (°C)\", color=MUTED_HEX, fontsize=11, labelpad=10)\n",
        "for spine in [\"top\", \"right\"]: ax.spines[spine].set_visible(False)\n",
        "ax.spines[\"left\"].set_color(GRID); ax.spines[\"bottom\"].set_color(GRID)\n",
        "ax.grid(axis=\"x\", linestyle=\"-\", color=GRID, alpha=0.8)\n",
        "ax.grid(axis=\"y\", visible=False)\n",
        "\n",
        "for bar in bars:\n",
        "    width = bar.get_width()\n",
        "    ax.text(width + 0.05, bar.get_y() + bar.get_height()/2, f\"+{width:.2f}\", va=\"center\", color=TEXT, fontsize=10, fontweight=\"bold\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "defceaaa"
      },
      "source": [
        "### Correlation: Heat vs. Rainfall\n",
        "A natural follow-up is whether warmer years also tend to be wetter or drier. Plotting the national annual mean temperature against total precipitation gives a first look at how the two aggregates move together."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 537
        },
        "id": "53ddb414",
        "outputId": "d31348a4-bef7-49dd-d362-d8a8f695e3b9"
      },
      "outputs": [],
      "source": [
        "import seaborn as sns\n",
        "\n",
        "# Prepare data for correlation\n",
        "corr_df = italy_temperature_trend[['Year', 'Temperature']].merge(\n",
        "    italy_precipitation_trend[['Year', 'Precipitation']],\n",
        "    on='Year'\n",
        ")\n",
        "\n",
        "fig, ax = plt.subplots(figsize=(8, 6))\n",
        "sns.regplot(data=corr_df, x='Temperature', y='Precipitation',\n",
        "            scatter_kws={'color': ACCENT, 's': 60},\n",
        "            line_kws={'color': TEXT, 'linestyle': '--', 'linewidth': 1.5}, ax=ax)\n",
        "\n",
        "ax.set_title(\"Correlation: Annual Temp vs. Precipitation\", loc='left', pad=25)\n",
        "ax.set_xlabel(\"Mean Temperature (°C)\")\n",
        "ax.set_ylabel(\"Total Precipitation (mm)\")\n",
        "\n",
        "# Calculate correlation coefficient\n",
        "r_val = corr_df['Temperature'].corr(corr_df['Precipitation'])\n",
        "ax.text(0.05, 0.05, f\"Pearson r = {r_val:.2f}\", transform=ax.transAxes,\n",
        "        fontsize=12, fontweight='bold', color=TEXT)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "200f0b4a"
      },
      "source": [
        "The scatter shows a weak-to-moderate negative association in these annual aggregates: across the 18 years, warmer years have tended to coincide with somewhat lower national precipitation. With a Pearson r of about −0.34 the relationship is real but far from tight — plenty of warm years were not notably dry. This describes only how the two national means co-vary from year to year; it does not say anything about the physical mechanism driving either series, and a single year can move the relationship substantially — the sharp 2022 dip, for instance, reflects the severe drought Italy faced that year and visibly pulls the regression line down.\n",
        "\n",
        "# Mapping the change\n",
        "\n",
        "A time-series chart shows *when* temperatures rise. A map shows *where* the warmth sits. To geocode the provincial capitals once, I cache the coordinates locally and reuse them on subsequent renders."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 717
        },
        "id": "9bb2c621",
        "outputId": "8e4d31dd-dbe4-426f-a3ab-6ea6ac492f88"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import plotly.express as px\n",
        "import plotly.io as pio\n",
        "\n",
        "# Re-defining the parser to ensure data is available\n",
        "def parse_table(sheet_name, value_name):\n",
        "    df = pd.read_excel(FILE_PATH, sheet_name=sheet_name, header=None)\n",
        "    years = [int(x) for x in df.iloc[3, 1:19].tolist()]\n",
        "    data = df.iloc[4:, :19].copy()\n",
        "    data.columns = [\"Province\"] + years\n",
        "    data = data[data[\"Province\"].notna()]\n",
        "    data[\"Province\"] = data[\"Province\"].astype(str).str.strip()\n",
        "    data = data[~data[\"Province\"].str.contains(r\"Fonte|Totale|Nota|Nelle tavole|^$\", regex=True, na=False)]\n",
        "    long = data.melt(id_vars=\"Province\", var_name=\"Year\", value_name=value_name)\n",
        "    long[\"Year\"] = pd.to_numeric(long[\"Year\"], errors=\"coerce\")\n",
        "    long[value_name] = pd.to_numeric(long[value_name].replace(\"....\", np.nan), errors=\"coerce\")\n",
        "    return long\n",
        "\n",
        "# Load data\n",
        "temp_df = parse_table(\"Tavola_1\", \"Temperature\")\n",
        "precip_df = parse_table(\"Tavola_2\", \"Precipitation\")\n",
        "italy_weather_data = temp_df.merge(precip_df, on=[\"Province\", \"Year\"], how=\"inner\")\n",
        "\n",
        "# A dictionary of coordinates for Italian provincial capitals\n",
        "city_coords = {\n",
        "    'Torino': [45.0703, 7.6869], 'Vercelli': [45.3238, 8.4232], 'Novara': [45.4468, 8.6212],\n",
        "    'Cuneo': [44.3833, 7.5500], 'Asti': [44.9005, 8.2069], 'Alessandria': [44.9129, 8.6154],\n",
        "    'Biella': [45.5630, 8.0579], 'Verbano-Cusio-Ossola': [45.9220, 8.5516], 'Aosta': [45.7371, 7.3201],\n",
        "    'Milano': [45.4642, 9.1900], 'Lodi': [45.3139, 9.5032], 'Monza e della Brianza': [45.5845, 9.2744],\n",
        "    'Bergamo': [45.6983, 9.6773], 'Brescia': [45.5398, 10.2181], 'Pavia': [45.1850, 9.1546],\n",
        "    'Como': [45.8081, 9.0852], 'Varese': [45.8167, 8.8333], 'Sondrio': [46.1690, 9.8731],\n",
        "    'Lecco': [45.8559, 9.3977], 'Mantova': [45.1564, 10.7911], 'Cremona': [45.1333, 10.0333],\n",
        "    'Bolzano': [46.4983, 11.3548], 'Trento': [46.0667, 11.1167], 'Verona': [45.4384, 10.9916],\n",
        "    'Vicenza': [45.5479, 11.5446], 'Belluno': [46.1408, 12.2161], 'Treviso': [45.6667, 12.2450],\n",
        "    'Venezia': [45.4408, 12.3155], 'Padova': [45.4064, 11.8768], 'Rovigo': [45.0711, 11.7905],\n",
        "    'Pordenone': [45.9569, 12.6605], 'Udine': [46.0625, 13.2346], 'Gorizia': [45.9409, 13.6222],\n",
        "    'Trieste': [45.6495, 13.7768], 'Piacenza': [45.0526, 9.6930], 'Parma': [44.8015, 10.3279],\n",
        "    'Reggio nell\\'Emilia': [44.6982, 10.6312], 'Modena': [44.6471, 10.9252], 'Bologna': [44.4949, 11.3426],\n",
        "    'Ferrara': [44.8381, 11.6198], 'Ravenna': [44.4183, 12.2035], 'Forlì-Cesena': [44.2227, 12.0409],\n",
        "    'Rimini': [44.0594, 12.5684], 'Massa Carrara': [44.0375, 10.1417], 'Lucca': [43.8429, 10.5027],\n",
        "    'Pistoia': [43.9333, 10.9167], 'Firenze': [43.7696, 11.2558], 'Livorno': [43.5485, 10.3106],\n",
        "    'Pisa': [43.7085, 10.4036], 'Arezzo': [43.4631, 11.8781], 'Siena': [43.3186, 11.3306],\n",
        "    'Grosseto': [42.7667, 11.1167], 'Prato': [43.8777, 11.1022], 'Perugia': [43.1107, 12.3908],\n",
        "    'Terni': [42.5639, 12.6427], 'Ancona': [43.6158, 13.5189], 'Pesaro e Urbino': [43.9100, 12.9133],\n",
        "    'Macerata': [43.3003, 13.4531], 'Ascoli Piceno': [42.8536, 13.5768], 'Fermo': [43.1610, 13.7183],\n",
        "    'Viterbo': [42.4191, 12.1051], 'Rieti': [42.4000, 12.8667], 'Roma': [41.9028, 12.4964],\n",
        "    'Latina': [41.4676, 12.9036], 'Frosinone': [41.6394, 13.3411], 'L\\'Aquila': [42.3489, 13.3980],\n",
        "    'Teramo': [42.6589, 13.7044], 'Pescara': [42.4618, 14.2142], 'Chieti': [42.3511, 14.1675],\n",
        "    'Campobasso': [41.5603, 14.6584], 'Isernia': [41.5947, 14.2342], 'Caserta': [41.0736, 14.3347],\n",
        "    'Benevento': [41.1297, 14.7821], 'Napoli': [40.8518, 14.2681], 'Avellino': [40.9144, 14.7936],\n",
        "    'Salerno': [40.6779, 14.7658], 'Foggia': [41.4622, 15.5446], 'Bari': [41.1171, 16.8719],\n",
        "    'Taranto': [40.4677, 17.2433], 'Brindisi': [40.6321, 17.9361], 'Lecce': [40.3515, 18.1750],\n",
        "    'Barletta': [41.3197, 16.2768], 'Potenza': [40.6404, 15.8051], 'Matera': [40.6664, 16.6043],\n",
        "    'Cosenza': [39.2983, 16.2537], 'Catanzaro': [38.9098, 16.5877], 'Reggio di Calabria': [38.1105, 15.6434],\n",
        "    'Crotone': [39.0808, 17.1273], 'Vibo Valentia': [38.6753, 16.1011], 'Messina': [38.1938, 15.5540],\n",
        "    'Palermo': [38.1157, 13.3615], 'Trapani': [38.0175, 12.5150], 'Agrigento': [37.3111, 13.5765],\n",
        "    'Caltanissetta': [37.4903, 14.0622], 'Enna': [37.5671, 14.2750], 'Catania': [37.5079, 15.0830],\n",
        "    'Ragusa': [36.9269, 14.7231], 'Siracusa': [37.0755, 15.2866], 'Sassari': [40.7259, 8.5615],\n",
        "    'Nuoro': [40.3231, 9.3303], 'Cagliari': [39.2238, 9.1217], 'Oristano': [39.9056, 8.5911],\n",
        "    'Carbonia': [39.1671, 8.5222],\n",
        "    'La Spezia': [44.2384, 9.6912],\n",
        "    'Andria': [41.2276, 16.2955],\n",
        "    'Trani': [41.2782, 16.4186],\n",
        "}\n",
        "\n",
        "# The workbook uses short province-capital names (e.g. \"Verbania\", \"Monza\",\n",
        "# \"Pesaro Urbino\", \"Forlì\") that do not match the longer administrative keys\n",
        "# in city_coords above. Normalize before lookup so no point disappears\n",
        "# silently when the dictionary key uses a different spelling.\n",
        "COORD_NAME_ALIASES = {\n",
        "    'Verbania': 'Verbano-Cusio-Ossola',\n",
        "    'Monza': 'Monza e della Brianza',\n",
        "    'Pesaro Urbino': 'Pesaro e Urbino',\n",
        "    'Forlì': 'Forlì-Cesena',\n",
        "}\n",
        "\n",
        "def lookup_coords(province):\n",
        "    key = COORD_NAME_ALIASES.get(province, province)\n",
        "    return city_coords.get(key, [None, None])\n",
        "\n",
        "# Map coordinates and warn about anything still unmapped so a miss is loud,\n",
        "# not silent.\n",
        "italy_weather_data['Lat'] = italy_weather_data['Province'].map(lambda x: lookup_coords(x)[0])\n",
        "italy_weather_data['Lon'] = italy_weather_data['Province'].map(lambda x: lookup_coords(x)[1])\n",
        "unmapped = italy_weather_data.loc[italy_weather_data['Lat'].isna(), 'Province'].unique()\n",
        "if len(unmapped) > 0:\n",
        "    print(f\"Unmapped provinces (no coordinates): {list(unmapped)}\")\n",
        "map_df = italy_weather_data.dropna(subset=['Lat', 'Lon', 'Temperature']).copy()\n",
        "\n",
        "# Prepare Mapbox map\n",
        "fig = px.scatter_mapbox(\n",
        "    map_df,\n",
        "    lat=\"Lat\", lon=\"Lon\", color=\"Temperature\",\n",
        "    hover_name=\"Province\", animation_frame=\"Year\",\n",
        "    title=\"Evolution of Annual Temperatures in Italy (2006-2023)\",\n",
        "    color_continuous_scale=px.colors.sequential.YlOrRd,\n",
        "    range_color=[map_df['Temperature'].min(), map_df['Temperature'].max()],\n",
        "    mapbox_style=\"open-street-map\",\n",
        "    zoom=4.5, center=dict(lat=42.0, lon=12.5)\n",
        ")\n",
        "\n",
        "fig.update_layout(height=700, margin=dict(r=10, t=80, b=40, l=10), paper_bgcolor=globals().get('BG', '#f7f6f2'))\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LzB-x7Gf-z21"
      },
      "source": [
        "A companion view averages each capital across the whole 2006–2023 window, so the persistent north–south gradient — Sicily and the far south glowing warmest, the Alps coolest — is easy to read at a glance."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 717
        },
        "id": "B_YHcZDT9u1X",
        "outputId": "23702d8e-935d-44b2-bf45-2070cd70cd16"
      },
      "outputs": [],
      "source": [
        "# Calculate mean per city for the static map\n",
        "mean_per_city = (\n",
        "    map_df.groupby(\"Province\", as_index=False)\n",
        "    .agg(\n",
        "        Lat=(\"Lat\", \"first\"),\n",
        "        Lon=(\"Lon\", \"first\"),\n",
        "        Temperature=(\"Temperature\", \"mean\"),\n",
        "    )\n",
        ")\n",
        "\n",
        "# Create a static Mapbox plot to match the interactive look and feel\n",
        "fig_static = px.scatter_mapbox(\n",
        "    mean_per_city,\n",
        "    lat=\"Lat\",\n",
        "    lon=\"Lon\",\n",
        "    color=\"Temperature\",\n",
        "    hover_name=\"Province\",\n",
        "    size_max=15,\n",
        "    color_continuous_scale=px.colors.sequential.YlOrRd,\n",
        "    range_color=[map_df['Temperature'].min(), map_df['Temperature'].max()],\n",
        "    mapbox_style=\"open-street-map\",\n",
        "    zoom=4.5,\n",
        "    center=dict(lat=42.0, lon=12.5),\n",
        "    title=\"<b>Mean temperature by provincial capital (2006–2023)</b><br><span style='font-size:12px; color:#7a7974'>Average temperature across the full series from ISTAT Tavola_1</span>\"\n",
        ")\n",
        "\n",
        "fig_static.update_layout(\n",
        "    height=700,\n",
        "    margin=dict(r=10, t=80, b=40, l=10),\n",
        "    paper_bgcolor=BG,\n",
        "    font=dict(color=TEXT)\n",
        ")\n",
        "\n",
        "fig_static.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4SBowc8X5c18"
      },
      "source": [
        "# Takeaways\n",
        "\n",
        "The updated workbook makes the national pattern easy to see: Italian provincial capitals show a clear upward temperature trend over 2006–2023, while precipitation is much noisier and less uniform.\n",
        "\n",
        "A few broad conclusions stand out:\n",
        "\n",
        "- The average across provincial capitals points to steady warming across the period.\n",
        "- The warmest cities are concentrated in Sicily and the far south.\n",
        "- The \"fastest warming\" ranking depends strongly on data completeness, so filtering sparse series matters.\n",
        "- The 2023 anomaly table adds important historical context, showing that several cities were not only warm in absolute terms, but unusually warm relative to climatological baselines.\n",
        "\n",
        "The natural next extension would be to pull in the regional-capital sheets (`Tavola_9` and `Tavola_10`) and the extremes tables (`Tavola_7`, `Tavola_8`, `Tavola_11`–`Tavola_14`) for a second article focused on heat extremes, warm nights, and heavy-rain indicators."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ABNK5BWy9dgF"
      },
      "source": [
        "The data, the original notebook, and this article are all open source. If this analysis was useful to you, a star on GitHub is the simplest way to say so — and it helps others find it too.\n",
        "\n",
        "[⭐ Star the repo on GitHub](https://github.com/feddernico/temperature-italia)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "reZYxCLW9lRT"
      },
      "source": [
        "Prefer to run it yourself? Grab the Jupyter notebook behind this post, complete with the ISTAT workbook bundled alongside, and reproduce every chart locally.\n",
        "\n",
        "[⬇️ Download the notebook](/downloads/italian-temperatures-2006-2023.ipynb)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "M3-L_tn39bo4"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
