Global temperature changes from 12,000 weather stations (1850-2026)

[sselvakumaran]

Dataset: https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network-monthly

[sselvakumaran]

code:

#!/usr/bin/env python3
"""
NOTE: this is a simplified version provided by Claude Opus 5, removing specific
styling features (rending options, fonts) that decreases reproducibility
and adding explainatory comments. enjoy!

Every station with a long record is drawn as one faint line showing 
how much warmer or cooler it is than its own first years.
~14.8k stations, ~1M line segments.

Data:  NOAA GHCN-Monthly v4 (adjusted), auto-downloaded on first run (~44 MB).
Needs: xy + numpy only. Run: python warming_trajectories_simple.py
"""
import glob, io, os, tarfile, urllib.request
import numpy as np
import xy

# --- 1. download NOAA's monthly station temperatures (only the first time) ---
DATA = "data/ghcn"
if not glob.glob(f"{DATA}/**/*.qcf.dat", recursive=True):
    os.makedirs(DATA, exist_ok=True)
    url = "https://www.ncei.noaa.gov/pub/data/ghcn/v4/ghcnm.tavg.latest.qcf.tar.gz"
    print("downloading GHCN-M v4 ...")
    blob = urllib.request.urlopen(url, timeout=120).read()
    tarfile.open(fileobj=io.BytesIO(blob), mode="r:gz").extractall(DATA)
dat = glob.glob(f"{DATA}/**/*.qcf.dat", recursive=True)[0]

# --- 2. read the fixed-width file (each row = one station-year) ---
# Layout per row: chars 0-10 station id, 11-14 year, then 12 monthly temperatures
# as 8-char blocks from char 19 (value in hundredths of a degree C, -9999 = missing).
rows = np.frombuffer(open(dat, "rb").read(), np.uint8)
rows = rows[: rows.size // 116 * 116].reshape(-1, 116)[:, :115]  # 115 chars + newline

def to_int(cols):
    # Decode a block of fixed-width ASCII digit columns into signed ints, all at
    # once (no slow Python loop): digit value * place value, negated if a '-' (45).
    digits = np.where((cols >= 48) & (cols <= 57), cols - 48, 0).astype(np.int64)
    value = digits @ (10 ** np.arange(cols.shape[1] - 1, -1, -1)).astype(np.int64)
    return np.where((cols == 45).any(1), -value, value)

year = to_int(rows[:, 11:15])
_, station = np.unique(np.ascontiguousarray(rows[:, :11]).view("S11").ravel(),
                       return_inverse=True)                    # station id -> 0,1,2,...
months = np.stack([to_int(rows[:, 19 + m * 8:24 + m * 8]) for m in range(12)], 1)

# --- 3. one number per station-year: the annual mean (need >= 6 good months) ---
good = months != -9999
annual = np.where(good, months, 0).sum(1) / np.maximum(good.sum(1), 1) / 100.0  # deg C
use = (good.sum(1) >= 6) & (year >= 1850)
station, year, annual = station[use], year[use], annual[use]

# sort so every station's years sit together in time order
order = np.lexsort((year, station))
station, year, annual = station[order], year[order], annual[order]
starts = np.concatenate([[True], station[1:] != station[:-1]])  # first row of a station
ends = np.concatenate([station[:-1] != station[1:], [True]])    # last row of a station

# --- 4. smooth each station over a 5-year window, then subtract its first value ---
i = np.arange(station.size)
lo = np.maximum(np.maximum.accumulate(np.where(starts, i, -1)), i - 2)          # window start
hi = np.minimum(np.minimum.accumulate(np.where(ends, i, i.size)[::-1])[::-1], i + 2)  # window end
csum = np.concatenate([[0.0], np.cumsum(annual)])
smooth = (csum[hi + 1] - csum[lo]) / (hi - lo + 1)              # 5-year running mean
first_val = np.full(station.max() + 1, np.nan)
first_val[station[starts]] = smooth[starts]
delta = smooth - first_val[station]                            # deg C vs the station's start

# --- 5. keep stations with 40+ years, then link each year to the next as a segment ---
long_enough = np.zeros(station.max() + 1, bool)
long_enough[station[ends]] = year[ends] - year[starts] >= 40
net = np.full(station.max() + 1, np.nan)                       # per-station total change
net[station[ends]] = delta[ends]                               # used for the line color
link = (station[:-1] == station[1:]) & long_enough[station[:-1]]

# --- 6. plot ~1M line segments, colored blue (cooled) to red (warmed) ---
chart = xy.scatter_chart(
    xy.segments(year[:-1][link], delta[:-1][link], year[1:][link], delta[1:][link],
                color=net[station[:-1][link]], colormap="rdbu_r",
                domain=(-2.5, 2.5), width=0.7, opacity=0.20),
    xy.x_axis(label="Year", grid=False),
    # y spans -5..5 so the tallest lines spill past the -4..4 ticks; 0 sits centered
    xy.y_axis(label="temperature vs. the station's first years (deg C)",
              domain=(-5, 5), label_position="center", grid=False, line=False,
              tick_values=[-4, -2, 0, 2, 4]),
    xy.theme(background="#000000", plot_background="#000000",
             text_color="#c9c9c9", axis_color="#555555"),
    title="Warming, one line per weather station",
    width=1920, height=1200,
)
chart.to_png("warming.png", scale=2)   # native renderer: no browser / fonts needed
print("saved warming.png")

[sselvakumaran]

also seen here! https://www.reddit.com/r/dataisbeautiful/comments/1vc2p7e/oc_global_temperature_changes_from_12000_weather/