9/11/2026

It was hot in Annapolis this summary (~early July 2026). I'm interested in how hot it was compared to other periods.

Data

To answer this question we need data!

A good dataset to use is ERA5-land given it's long time series (1950 to present). 

The European Commission via Destination Earth hosts the data on the Earth Data Hub.

Query

You can query the data for a single point for Annapolis doing uv run etl.py

etl.py:

# /// script

# dependencies = [

#   "xarray",

#   "zarr",

#   "pandas",

#   "geopy",

#   "pyarrow",

#   "fsspec",

#   "aiohttp",

#   "dask",

#   "numpy",

# ]

# ///

import os

import numpy as np

import pandas as pd

import xarray as xr

from geopy.geocoders import Photon, Nominatim

from geopy.adapters import GeocoderUnavailable

# Connect to Earth Data Hub ERA5 Land daily dataset

api_key = os.environ.get("EARTHDATAHUB_API_KEY")

zarr_url = f"https://edh:{api_key}@api.earthdatahub.destine.eu/era5/era5-land-daily-utc-v1.zarr"

ds = xr.open_dataset(

    zarr_url,

    chunks={},

    engine="zarr",

    zarr_format=3,

)

# Geocode Annapolis coordinates

city = "Annapolis"

try:

    location = Photon().geocode(city)

except GeocoderUnavailable:

    location = Nominatim(user_agent="_").geocode(city)

lat, lon = location.point.latitude, location.point.longitude

lon_360 = lon + 360

# Search for the nearest valid land grid point (non-NaN)

# ERA5 Land resolution is 0.1 degrees. Search in a 0.3 degree radius.

search_range = 0.3

ds_search = ds.sel(

    latitude=slice(lat + search_range, lat - search_range),

    longitude=slice(lon_360 - search_range, lon_360 + search_range)

)

# Load a sample from the beginning of the time series to identify land grid cells (non-NaN)

sample = ds_search["t2m"].isel(valid_time=0).compute()

valid_mask = ~np.isnan(sample.values)

lats = sample.latitude.values

lons = sample.longitude.values

if not np.any(valid_mask):

    raise ValueError("No valid land grid cells found in the search area!")

best_lat, best_lon = None, None

min_dist = float("inf")

for r in range(len(lats)):

    for c in range(len(lons)):

        if valid_mask[r, c]:

            dist = (lats[r] - lat)**2 + (lons[c] - lon_360)**2

            if dist < min_dist:

                min_dist = dist

                best_lat = lats[r]

                best_lon = lons[c]

print(f"Nearest valid land cell found at: Lat {best_lat:.4f}, Lon {best_lon:.4f} (distance: {np.sqrt(min_dist):.4f} degrees)")

# Select the full time series for this location

ds_local = ds.sel(latitude=best_lat, longitude=best_lon)

# Extract temperature variable and load to pandas

ds_t2m = ds_local["t2m"]

df_t2m = ds_t2m.to_pandas()

# Convert Series to DataFrame if necessary, and write to a Parquet file

if isinstance(df_t2m, pd.Series):

    df_t2m = df_t2m.to_frame(name="t2m")

output_file = "era5_land_t2m_annapolis.parquet"

df_t2m.to_parquet(output_file)

print(f"Successfully saved ERA5 data to {output_file}")

 As of 9/7112026 the data goes from 1950-01-01 to 2026-08-31.

We can convert the daily time series into percentiles 

df_t2m = ds_t2m.to_pandas()

df_t2m["percentile"] = df_t2m["t2m"].rank(pct=True) * 100

 any find the max percentil this summer

df_t2m.reset_index()[df_t2m.reset_index()["valid_time"].dt.year == 2026]["percentile"].max()

 99.978!

That means there was only 6 hotter days than July 3rd 2026 since 1950 ( 28,002 days)

I knew it was a hot day!

# /// script

# dependencies = ["pandas", "pyarrow", "numpy", "matplotlib"]

# ///

import numpy as np, pandas as pd

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt

import matplotlib.dates as mdates

from matplotlib.patches import FancyBboxPatch

from matplotlib.lines import Line2D

SRC = "/mnt/c/Users/ray.bell/Downloads/annapolis_t2m_era5land.parquet"

C2F = lambda c: c * 9 / 5 + 32

df = pd.read_parquet(SRC)

df["C"] = df.t2m - 273.15

df["F"] = C2F(df.C)

df["year"] = df.valid_time.dt.year

df["md"] = df.valid_time.dt.strftime("%m-%d")

# ---- headline stats -------------------------------------------------

j3 = df[df.valid_time == "2026-07-03"].iloc[0]

n_hotter = int((df.C > j3.C).sum())

rank_j3 = n_hotter + 1

jja = df[df.valid_time.dt.month.isin([6, 7, 8])]

smean = jja.groupby("year")["F"].mean().sort_values(ascending=False)

rank_2026 = list(smean.index).index(2026) + 1

base = jja[jja.year <= 2020].F.mean()

anom = smean[2026] - base

# ---- day-of-year climatology, 1950-2025, +/-7 day window ------------

hist = df[df.year <= 2025].copy()

hist["doy"] = hist.valid_time.dt.dayofyear

# normalise leap years so Mar-Dec aligns

leap = hist.valid_time.dt.is_leap_year & (hist.doy > 59)

hist.loc[leap, "doy"] -= 1

by_doy = {d: g.F.values for d, g in hist.groupby("doy")}

s26 = df[(df.year == 2026) & df.valid_time.dt.month.isin([6, 7, 8])].copy()

s26["doy"] = s26.valid_time.dt.dayofyear  # 2026 not a leap year

rows = []

for d in s26.doy:

    w = np.concatenate([by_doy[((d + o - 1) % 365) + 1] for o in range(-7, 8)])

    rows.append([np.percentile(w, p) for p in (1, 10, 90, 99)])

clim = np.array(rows)

p01, p10, p90, p99 = clim[:, 0], clim[:, 1], clim[:, 2], clim[:, 3]

x = s26.valid_time.values

y = s26.F.values

print(f"Jul 3: {j3.F:.1f}F / {j3.C:.2f}C  rank {rank_j3} of {len(df):,}  ({n_hotter} hotter)")

print(f"JJA 2026 mean {smean[2026]:.1f}F  rank {rank_2026} of {len(smean)}  anom +{anom:.1f}F")

print(f"days above the 90th pct for the date: {int((y > p90).sum())} of {len(y)}")

print(f"days above the 99th pct for the date: {int((y > p99).sum())} of {len(y)}")

# ---- theme ----------------------------------------------------------

THEMES = {

    "light": dict(surface="#fcfcfb", page="#f9f9f7", ink="#0b0b0b", ink2="#52514e",

                  muted="#898781", grid="#e1e0d9", axis="#c3c2b7",

                  accent="#eb6834", band1="#dedcd4", band2="#efeee9", suffix=""),

    "dark":  dict(surface="#1a1a19", page="#0d0d0d", ink="#ffffff", ink2="#c3c2b7",

                  muted="#898781", grid="#2c2c2a", axis="#383835",

                  accent="#d95926", band1="#33332f", band2="#232321", suffix="_dark"),

}

SANS = ["Inter", "system-ui", "DejaVu Sans"]

def build(t):

    plt.rcParams.update({"font.family": "sans-serif", "font.sans-serif": SANS,

                         "svg.fonttype": "none"})

    fig = plt.figure(figsize=(12.0, 7.4), dpi=200, facecolor=t["page"])

    fig.patch.set_facecolor(t["page"])

    # card surface

    card = FancyBboxPatch((0.012, 0.012), 0.976, 0.976, transform=fig.transFigure,

                          boxstyle="round,pad=0,rounding_size=0.012",

                          facecolor=t["surface"], edgecolor="none", zorder=-10)

    fig.add_artist(card)

    L, R = 0.068, 0.941

    # ---------------- headline text ----------------

    fig.text(L, 0.945, "Annapolis just had one of its hottest summers since 1950",

             ha="left", va="top", fontsize=21, color=t["ink"], weight="semibold")

    fig.text(L, 0.898,

             "Daily mean air temperature, June–August 2026, against the range for "

             "the same date over 1950–2025.",

             ha="left", va="top", fontsize=11.5, color=t["ink2"])

    # ---------------- hero + stat row ----------------

    hy = 0.792

    fig.text(L, hy, f"{n_hotter}", ha="left", va="center", fontsize=62,

             color=t["ink"], weight="bold")

    hw = 0.052 if n_hotter < 10 else 0.088

    fig.text(L + hw, hy + 0.028, "days hotter than 3 July 2026",

             ha="left", va="center", fontsize=12.5, color=t["ink"], weight="semibold")

    fig.text(L + hw, hy - 0.016, f"in the {len(df):,} days since 1 January 1950",

             ha="left", va="center", fontsize=11, color=t["muted"])

    tiles = [(f"{j3.F:.1f}°F", "3 July daily mean", f"{j3.C:.1f}°C · rank {rank_j3} of {len(df):,}"),

             (f"{smean[2026]:.0f}°F", "Summer 2026 average", f"{rank_2026}th-hottest of {len(smean)} summers"),

             (f"+{anom:.1f}°F", "vs the 1950–2020 norm", "for June–August")]

    tx = 0.455

    for val, lab, sub in tiles:

        fig.text(tx, hy + 0.045, lab, ha="left", va="center", fontsize=10,

                 color=t["muted"])

        fig.text(tx, hy + 0.004, val, ha="left", va="center", fontsize=25,

                 color=t["ink"], weight="semibold")

        fig.text(tx, hy - 0.042, sub, ha="left", va="center", fontsize=9.5,

                 color=t["ink2"])

        tx += 0.178

    fig.add_artist(Line2D([L, R], [0.706, 0.706], color=t["grid"], lw=1,

                          transform=fig.transFigure))

    # ---------------- plot ----------------

    ax = fig.add_axes([L, 0.115, R - L, 0.522])

    ax.set_facecolor(t["surface"])

    for s in ax.spines.values():

        s.set_visible(False)

    ax.spines["bottom"].set_visible(True)

    ax.spines["bottom"].set_color(t["axis"])

    ax.spines["bottom"].set_linewidth(1)

    ax.fill_between(x, p01, p99, color=t["band2"], lw=0, zorder=1)

    ax.fill_between(x, p10, p90, color=t["band1"], lw=0, zorder=2)

    ax.set_axisbelow(True)

    ax.grid(axis="y", color=t["grid"], lw=1, ls="-", zorder=0)

    ax.tick_params(axis="both", length=0, colors=t["muted"], labelsize=10.5)

    ax.plot(x, y, color=t["accent"], lw=2, solid_joinstyle="round",

            solid_capstyle="round", zorder=4)

    # July 3 marker: dot + 2px surface ring

    ax.plot([j3.valid_time], [j3.F], "o", ms=9, mfc=t["accent"],

            mec=t["surface"], mew=2, zorder=6)

    lo = min(p01.min(), y.min()) - 1.0

    hi = max(p99.max(), y.max()) + 5.0

    ax.set_ylim(lo, hi)

    ax.set_xlim(x[0], x[-1])

    # direct label on the extreme, with a leader line

    ax.annotate(f"3 July  {j3.F:.1f}°F",

                xy=(j3.valid_time, j3.F), xytext=(-14, 26),

                textcoords="offset points", ha="right", va="bottom",

                fontsize=11.5, color=t["ink"], weight="semibold",

                arrowprops=dict(arrowstyle="-", color=t["muted"], lw=1,

                                shrinkA=0, shrinkB=6))

    ax.xaxis.set_major_locator(mdates.MonthLocator())

    ax.xaxis.set_major_formatter(mdates.DateFormatter("%B"))

    ax.xaxis.set_minor_locator(mdates.DayLocator(bymonthday=(8, 15, 22, 29)))

    for lab in ax.get_xticklabels():

        lab.set_ha("left")

    ax.set_yticks(np.arange(60, 96, 5))

    ax.set_yticklabels([f"{v}°F" for v in np.arange(60, 96, 5)])

    # secondary Celsius axis (same measure, second unit - not a second scale)

    ax2 = ax.twinx()

    ax2.set_ylim(lo, hi)

    for s in ax2.spines.values():

        s.set_visible(False)

    cticks = np.arange(15, 34, 5)

    ax2.set_yticks(C2F(cticks))

    ax2.set_yticklabels([f"{v}°C" for v in cticks])

    ax2.tick_params(axis="y", length=0, colors=t["muted"], labelsize=9.5)

    leg = ax.legend(handles=[

        Line2D([], [], color=t["accent"], lw=2, label="2026"),

        Line2D([], [], color=t["band1"], lw=9, label="10th–90th percentile for the date"),

        Line2D([], [], color=t["band2"], lw=9, label="1st–99th percentile"),

    ], loc="lower left", bbox_to_anchor=(0.0, 1.005), frameon=False,

        fontsize=10.5, labelcolor=t["ink2"], handlelength=1.6,

        handletextpad=0.7, borderaxespad=0.0, ncol=3, columnspacing=2.0)

    for txt in leg.get_texts():

        txt.set_color(t["ink2"])

    fig.text(L, 0.045,

             "Data: ERA5-Land 2 m air temperature, daily mean, Annapolis, Maryland · "

             "1 Jan 1950 – 31 Aug 2026 (28,002 days)",

             ha="left", va="center", fontsize=9, color=t["muted"])

    out = f"/home/raybell/tmp/annapolis_summer_2026{t['suffix']}.png"

    fig.savefig(out, dpi=200, facecolor=t["page"])

    plt.close(fig)

    print("wrote", out)

for t in THEMES.values():

    build(t)

# ---- data twin: every plotted value in text form ---------------------

tw = pd.DataFrame({

    "date": s26.valid_time.dt.strftime("%Y-%m-%d"),

    "temp_F": np.round(y, 1), "temp_C": np.round(s26.C.values, 2),

    "clim_p10_F": np.round(p10, 1), "clim_p90_F": np.round(p90, 1),

    "clim_p01_F": np.round(p01, 1), "clim_p99_F": np.round(p99, 1),

})

tw["above_p90_for_date"] = tw.temp_F > tw.clim_p90_F

tw.to_csv("/home/raybell/tmp/annapolis_summer_2026_data.csv", index=False)

print("wrote /home/raybell/tmp/annapolis_summer_2026_data.csv")