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")