How have strong El Niño winters changed U.S. snowfall, and what does that suggest for 2026–27?

from datetime import UTC, datetime
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import usdata
from usdata import cite_lockfile, pull, verify

# One figure style for every usdata notebook, so previews look alike.
plt.rcParams.update(
    {
        "figure.figsize": (8, 4.5),
        "figure.dpi": 120,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.grid": True,
        "grid.alpha": 0.3,
        "font.size": 10,
    }
)
manifest = Path("dataset.yaml")
print("Executed (UTC):", datetime.now(UTC).isoformat(timespec="seconds"))
print(f"usdata {usdata.__version__}; pandas {pd.__version__}")

REGIONS = {  # west, south, east, north in degrees; each box is half-open
    "pacific-northwest": (-125.0, 42.0, -116.5, 49.0),
    "northern-rockies": (-116.5, 43.0, -104.0, 49.0),
    "upper-midwest": (-97.5, 43.5, -87.0, 49.0),
    "great-lakes-ohio-valley": (-89.0, 38.0, -78.0, 43.5),
    "northeast": (-78.0, 36.5, -67.0, 47.5),
    "southern-plains": (-104.0, 32.0, -94.0, 37.0),
    "southwest": (-114.0, 31.3, -104.0, 37.0),
}
result = pull(manifest)
for region in REGIONS:
    items = result.by_source[region]
    size = sum(item.provenance.size for item in items)
    print(f"{region:24} {len(items):3} stations  {size / 1e6:5.1f} MB")
print("Retrieved (UTC):", result.one("roni").provenance.retrieved_at.isoformat(timespec="seconds"))

FIRST, LAST, STRONG = 1949, 2025, 1.5


def winter_peaks(name: str) -> pd.Series:
    """The peak OND, NDJ, or DJF value of one CPC index, per winter."""
    table = pd.read_csv(result.one(name).path, sep=r"\s+")
    table = table[table["SEAS"].isin(["OND", "NDJ", "DJF"])]
    winter = table["YR"] - (table["SEAS"] == "DJF")
    return table["ANOM"].groupby(winter).max().loc[FIRST:LAST]


peaks = pd.DataFrame({"RONI": winter_peaks("roni"), "ONI": winter_peaks("oni")})
strong = {name: sorted(peaks.index[peaks[name] >= STRONG]) for name in peaks}
table = peaks.loc[sorted(set(strong["RONI"]) | set(strong["ONI"]))]
table.index = [f"{year}-{str(year + 1)[2:]}" for year in table.index]
table

frames = []
for region in REGIONS:
    for item in result.by_source[region]:
        frame = pd.read_csv(
            item.path,
            usecols=["STATION", "DATE", "LATITUDE", "LONGITUDE", "ELEVATION", "SNOW"],
            dtype={"STATION": "string"},
        )
        frames.append(frame.assign(region=region))
monthly = pd.concat(frames, ignore_index=True)
month = monthly["DATE"].str[5:].astype(int)
monthly = monthly[month.isin([11, 12, 1, 2, 3])].copy()
monthly["winter"] = monthly["DATE"].str[:4].astype(int) - (monthly["DATE"].str[5:].astype(int) <= 3)
monthly = monthly[monthly["winter"].between(FIRST, LAST)]
stations = monthly.groupby("STATION")[["region", "LATITUDE", "LONGITUDE", "ELEVATION"]].first()

counts = monthly.groupby(["STATION", "winter"])["SNOW"].agg(["count", "sum"])
winters = counts["sum"].where(counts["count"] == 5).unstack("winter") / 10  # cm
complete = winters.notna().sum(axis=1)
complete_strong = winters[strong["RONI"]].notna().sum(axis=1)
usual = winters.median(axis=1)
keep = (complete >= 60) & (complete_strong >= 6) & (usual >= 10)
long_record = complete >= 60
enough_strong = long_record & (complete_strong >= 6)
print(f"{len(winters)} stations; kept {int(keep.sum())}")
print(f"  fewer than 60 complete winters: {int((~long_record).sum())}")
print(f"  fewer than 6 complete strong winters: {int((long_record & ~enough_strong).sum())}")
print(f"  median winter under 10 cm: {int((enough_strong & (usual < 10)).sum())}")
snow = winters[keep]
stations = stations.loc[snow.index]
stations.groupby("region").size().rename("stations").to_frame().T

def ratios(table: pd.DataFrame, chosen) -> pd.Series:
    """Each station's median snowfall in the chosen winters over its median in all."""
    return table[list(chosen)].median(axis=1) / table.median(axis=1)


stations["ratio"] = ratios(snow, strong["RONI"])
regions = stations.groupby("region")["ratio"].agg(
    stations="size", median_ratio="median", share_below=lambda r: (r < 1).mean()
)
rng = np.random.default_rng(1950)
draws = np.array(
    [
        ratios(snow, rng.choice(snow.columns, size=len(strong["RONI"]), replace=False))
        .groupby(stations["region"])
        .median()
        .reindex(regions.index)
        for _ in range(2000)
    ]
)
distance = np.abs(np.log(regions["median_ratio"].to_numpy()))
regions["p"] = (np.abs(np.log(draws)) >= distance).mean(axis=0)
by_winter = (
    snow[strong["RONI"]].div(snow.median(axis=1), axis=0).groupby(stations["region"]).median()
)
regions["strong_winters_below"] = (by_winter < 1).sum(axis=1).astype(int)
regions = regions.loc[list(REGIONS)]
regions.round(2)

fig, ax = plt.subplots(layout="constrained")
mesh = ax.scatter(
    stations["LONGITUDE"],
    stations["LATITUDE"],
    c=np.log2(stations["ratio"]),
    cmap="BrBG",
    vmin=-1.5,
    vmax=1.5,
    s=22,
    edgecolor="0.3",
    linewidth=0.3,
)
for west, south, east, north in REGIONS.values():
    ax.add_patch(
        plt.Rectangle(
            (west, south),
            east - west,
            north - south,
            fill=False,
            edgecolor="0.6",
            linewidth=0.6,
            linestyle="--",
        )
    )
ax.set_aspect(1 / np.cos(np.radians(40)))
ax.grid(False)
ax.set_xlabel("Longitude (degrees east)")
ax.set_ylabel("Latitude (degrees north)")
ax.set_title("Snowfall in strong El Niño winters, relative to each station's usual")
bar = fig.colorbar(mesh, ax=ax, shrink=0.75, ticks=np.log2([0.35, 0.5, 0.71, 1, 1.41, 2, 2.83]))
bar.ax.set_yticklabels(["35%", "50%", "71%", "100%", "141%", "200%", "283%"])
bar.set_label("Median in the 8 strong winters / median in all 77")
plt.show()

last_winter = (snow[2023] / snow.median(axis=1)).groupby(stations["region"]).median()
comparison = pd.DataFrame(
    {
        "strong El Niño winters": regions["median_ratio"],
        "2023-24": last_winter.reindex(regions.index),
    }
)
comparison.round(2)

with_2023 = stations.assign(ratio_oni=ratios(snow, strong["ONI"]))
with_2023.groupby("region")[["ratio", "ratio_oni"]].median().loc[list(REGIONS)].round(2)

roni = pd.read_csv(result.one("roni").path, sep=r"\s+")
oni = pd.read_csv(result.one("oni").path, sep=r"\s+")
latest = roni.iloc[-1]
print(
    f"Latest season: {latest['SEAS']} {latest['YR']}: RONI {latest['ANOM']:+.2f}, "
    f"ONI {oni.iloc[-1]['ANOM']:+.2f}"
)
print("\nIf 2026-27 is a strong El Niño winter, the eight past ones suggest:")
for region, row in regions.iterrows():
    print(
        f"  {region:24} median {row['median_ratio']:.0%} of usual; "
        f"{int(row['strong_winters_below'])} of 8 strong winters below usual"
    )

assert verify(manifest) == []
for citation in cite_lockfile(manifest):
    print(citation.as_text())