Heat waves & public health — apparent temperature¶

Heat-related health risk depends on more than air temperature: humidity amplifies the physiological stress because evaporative cooling becomes less effective. The heat index (apparent temperature) combines air temperature and dewpoint into a single number that maps to discomfort categories.

Domain context. A common formula (Rothfusz, NOAA NWS) is

$$ HI = -42.379 + 2.04901523 T + 10.14333127 RH - 0.22475541 \, T \cdot RH - \dots $$

where $T$ is air temperature (°F) and $RH$ is relative humidity (%). Below $T \approx 80$°F the formula returns roughly $T$ itself; above it, humidity drives the index up sharply. A simpler proxy used here: the wet-bulb temperature, derivable from $T_{2m}$ and dewpoint $T_d$ via Davies-Jones / Stull approximations.

Setup¶

All the imports in one place. earthlens provides the unified EarthLens entry point plus AggregationConfig (to fold hourly CDS data down to daily means); pyramids supplies Dataset for reading the written GeoTIFFs; numpy / pandas / matplotlib do the maths and the plot.

In [ ]:

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pyramids.dataset import Dataset

from earthlens import AggregationConfig, EarthLens

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.dataset import Dataset from earthlens import AggregationConfig, EarthLens

Step 1 — pull a summer's daily data¶

Box: 1° around Athens, Greece (37°–38°N, 23°–24°E) — well-known European heat-wave hotspot. Range June–August 2022 covers the infamous summer heat wave. Daily resolution captures the day-to-day variability that monthly aggregation would smooth out.

First create the output directory and build the ECMWF request — the two ERA5 single-level variables (2m-temperature and 2m-dewpoint-temperature) over the Athens box for the summer window. Building the request is kept separate from running it.

In [ ]:

OUT = Path("data/era5-athens-summer")
OUT.mkdir(parents=True, exist_ok=True)

earthlens = EarthLens(
    data_source="ecmwf",
    cadence="daily",
    start="2022-06-01",
    end="2022-08-31",
    dataset="reanalysis-era5-single-levels",
    variables=[
        "2m-temperature",
        "2m-dewpoint-temperature",
    ],
    aoi=[23.0, 37.0, 24.0, 38.0],
    path=str(OUT),
)

OUT = Path("data/era5-athens-summer") OUT.mkdir(parents=True, exist_ok=True) earthlens = EarthLens( data_source="ecmwf", cadence="daily", start="2022-06-01", end="2022-08-31", dataset="reanalysis-era5-single-levels", variables=[ "2m-temperature", "2m-dewpoint-temperature", ], aoi=[23.0, 37.0, 24.0, 38.0], path=str(OUT), )

Run the download. Passing aggregate= chains the aggregator into the same call: the hourly CDS NetCDFs are folded to daily means (freq="1D", op="mean") and written under OUT/aggregated.

In [ ]:

earthlens.download(aggregate=AggregationConfig(freq="1D", op="mean"))

earthlens.download(aggregate=AggregationConfig(freq="1D", op="mean"))

Step 2 — derive RH and wet-bulb from T and Td¶

Magnus-Tetens for saturation vapour pressure, then RH = e/es. Stull (2011) for wet-bulb temperature.

Read the aggregated rasters¶

Stack the per-day daily-mean GeoTIFFs for each variable into a (day, y, x) array and convert from kelvin to °C.

In [ ]:

agg = OUT / "aggregated"
t_air = (
    np.stack(
        [
            Dataset.read_file(str(p)).read_array()
            for p in sorted(agg.glob("2m_temperature_1D_*.tif"))
        ]
    )
    - 273.15
)
t_dew = (
    np.stack(
        [
            Dataset.read_file(str(p)).read_array()
            for p in sorted(agg.glob("2m_dewpoint_temperature_1D_*.tif"))
        ]
    )
    - 273.15
)

agg = OUT / "aggregated" t_air = ( np.stack( [ Dataset.read_file(str(p)).read_array() for p in sorted(agg.glob("2m_temperature_1D_*.tif")) ] ) - 273.15 ) t_dew = ( np.stack( [ Dataset.read_file(str(p)).read_array() for p in sorted(agg.glob("2m_dewpoint_temperature_1D_*.tif")) ] ) - 273.15 )

Relative humidity from the Magnus formula¶

Saturation vapour pressure es(t) (hPa) from the Magnus formula; relative humidity is then the ratio of the dewpoint and air-temperature saturation pressures.

In [ ]:

# Magnus formula for saturation vapour pressure (hPa).


def es(t_celsius):
    return 6.112 * np.exp(17.67 * t_celsius / (t_celsius + 243.5))


rh = 100.0 * es(t_dew) / es(t_air)

# Magnus formula for saturation vapour pressure (hPa). def es(t_celsius): return 6.112 * np.exp(17.67 * t_celsius / (t_celsius + 243.5)) rh = 100.0 * es(t_dew) / es(t_air)

Wet-bulb temperature (Stull 2011)¶

The Stull (2011) closed-form approximation gives wet-bulb temperature from air temperature and RH; it is valid for rh > 5% and t in -20..50 °C.

In [ ]:

# Stull (2011) wet-bulb approximation, valid for rh > 5% and t -20..50 °C.


def wet_bulb(t, rh):
    return (
        t * np.arctan(0.151977 * (rh + 8.313659) ** 0.5)
        + np.arctan(t + rh)
        - np.arctan(rh - 1.676331)
        + 0.00391838 * (rh**1.5) * np.arctan(0.023101 * rh)
        - 4.686035
    )


t_wet_bulb = wet_bulb(t_air, rh)

# Stull (2011) wet-bulb approximation, valid for rh > 5% and t -20..50 °C. def wet_bulb(t, rh): return ( t * np.arctan(0.151977 * (rh + 8.313659) ** 0.5) + np.arctan(t + rh) - np.arctan(rh - 1.676331) + 0.00391838 * (rh**1.5) * np.arctan(0.023101 * rh) - 4.686035 ) t_wet_bulb = wet_bulb(t_air, rh)

Daily site-mean summary¶

Collapse each field over the spatial box to one value per day and assemble a tidy DataFrame indexed by date.

In [ ]:

site_t_air = np.nanmean(t_air, axis=(1, 2))
site_rh = np.nanmean(rh, axis=(1, 2))
site_t_wet_bulb = np.nanmean(t_wet_bulb, axis=(1, 2))

days = pd.date_range("2022-06-01", periods=len(site_t_air), freq="D")
summary = pd.DataFrame(
    {
        "T_2m [°C]": site_t_air.round(1),
        "RH [%]": site_rh.round(0),
        "T_wet-bulb [°C]": site_t_wet_bulb.round(1),
    },
    index=days,
)
summary.head()

site_t_air = np.nanmean(t_air, axis=(1, 2)) site_rh = np.nanmean(rh, axis=(1, 2)) site_t_wet_bulb = np.nanmean(t_wet_bulb, axis=(1, 2)) days = pd.date_range("2022-06-01", periods=len(site_t_air), freq="D") summary = pd.DataFrame( { "T_2m [°C]": site_t_air.round(1), "RH [%]": site_rh.round(0), "T_wet-bulb [°C]": site_t_wet_bulb.round(1), }, index=days, ) summary.head()

Step 3 — plot the heat-stress time series¶

Wet-bulb temperatures > 30 °C are physiologically stressful for healthy adults; > 35 °C is the theoretical limit of human survival (no evaporative cooling possible at all).

Plot the daily air and wet-bulb temperatures together, with the 30 °C and 35 °C wet-bulb thresholds drawn in and the stressful band shaded.

In [ ]:

fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(days, site_t_air, label="Air T_2m", color="tab:red", lw=1.5)
ax.plot(days, site_t_wet_bulb, label="Wet-bulb T", color="tab:purple", lw=2)
ax.axhline(30, color="orange", lw=0.8, ls="--", label="30 °C wet-bulb (stressful)")
ax.axhline(35, color="red", lw=0.8, ls="--", label="35 °C wet-bulb (survival limit)")
ax.fill_between(
    days,
    30,
    np.maximum(site_t_wet_bulb, 30),
    where=site_t_wet_bulb >= 30,
    color="orange",
    alpha=0.25,
)
ax.set_ylabel("Temperature [°C]")
ax.set_title("Athens summer 2022 — daily-mean air and wet-bulb temperatures")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(11, 5)) ax.plot(days, site_t_air, label="Air T_2m", color="tab:red", lw=1.5) ax.plot(days, site_t_wet_bulb, label="Wet-bulb T", color="tab:purple", lw=2) ax.axhline(30, color="orange", lw=0.8, ls="--", label="30 °C wet-bulb (stressful)") ax.axhline(35, color="red", lw=0.8, ls="--", label="35 °C wet-bulb (survival limit)") ax.fill_between( days, 30, np.maximum(site_t_wet_bulb, 30), where=site_t_wet_bulb >= 30, color="orange", alpha=0.25, ) ax.set_ylabel("Temperature [°C]") ax.set_title("Athens summer 2022 — daily-mean air and wet-bulb temperatures") ax.legend() ax.grid(alpha=0.3) plt.tight_layout() plt.show()

Notes¶

• Daily mean understates peak heat. Real heat-wave health exposure depends on daily maximum wet-bulb. To get it, retrieve hourly data (temporal_resolution="daily") and run aggregate_netcdf with op="max" instead of op="mean".
• CDS-Beta has hourly products. Use the derived-era5-single- levels-daily-statistics family for pre-aggregated daily max / min statistics — but those NetCDFs are already daily aggregates, so pass op="mean" explicitly when re-aggregating to coarser steps to avoid double-summing.
• Health categories. US NWS heat index categories: Caution > 27 °C, Extreme Caution > 32 °C, Danger > 41 °C, Extreme Danger > 54 °C (apparent temperature, not wet-bulb). The Stull wet-bulb here is an alternative — closer to the thermodynamic limit on evaporative cooling.
• Stull approximation limits. Valid for $T \in [-20, 50]$°C and $RH > 5\%$. Outside that range use the iterative Davies-Jones (2008) formulation.