|
#!/usr/bin/env -S uv run --script |
|
# |
|
# /// script |
|
# requires-python = ">=3.12" |
|
# dependencies = [ |
|
# "numpy>=2.5.1", |
|
# "typer>=0.27.0", |
|
# "xy>=0.0.5", |
|
# ] |
|
# /// |
|
|
|
# Original source: https://github.com/reflex-dev/xy/discussions/419 |
|
|
|
""" |
|
Fork of a script for plotting temperature changes from fist measurement |
|
|
|
NOTE: Original source was partially AI generated. |
|
|
|
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: `uv` only (pulling in xy, numpy, typer and their dependencies). |
|
|
|
Run `uv run tempvis.py` (or `uv run tempvis.py --help` for options). |
|
""" |
|
|
|
import glob |
|
import io |
|
import os |
|
import tarfile |
|
import urllib.request |
|
|
|
import numpy as np |
|
import xy |
|
|
|
# Defaults: |
|
DATA_DIR = "data/ghcn" |
|
TITLE = "Warming, one line per weather station" |
|
X_LABEL = "year" |
|
Y_LABEL = "temperature vs. the station's first years (deg C)" |
|
X_GRID = True |
|
Y_GRID = True |
|
URL = "https://www.ncei.noaa.gov/pub/data/ghcn/v4/ghcnm.tavg.latest.qcf.tar.gz" |
|
|
|
|
|
def get_rows( |
|
data: str = DATA_DIR, |
|
url: str = URL, |
|
force_download: bool = False, |
|
): |
|
data = "data/ghcn" |
|
globstr = f"{data}/**/*.qcf.dat" |
|
if not force_download and (existing := glob.glob(globstr, recursive=True)): |
|
path = existing[0] |
|
else: |
|
# --- 1. download NOAA's monthly station temperatures (only the first time) --- |
|
os.makedirs(data, exist_ok=True) |
|
print("downloading GHCN-M v4 ...", end="") |
|
blob = urllib.request.urlopen(url, timeout=120).read() |
|
with tarfile.open(fileobj=io.BytesIO(blob), mode="r:gz") as f: |
|
f.extractall(data) |
|
path = glob.glob(globstr, recursive=True)[0] |
|
print(" done.") |
|
|
|
with open(path, "rb") as f: |
|
rows = np.frombuffer(f.read(), np.uint8) |
|
|
|
# --- 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). |
|
return 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) |
|
|
|
|
|
def get_chart( |
|
title: str = TITLE, |
|
x_label: str = X_LABEL, |
|
y_label: str = Y_LABEL, |
|
show_x_grid: bool = X_GRID, |
|
show_y_grid: bool = Y_GRID, |
|
url: str = URL, |
|
data_dir: str = DATA_DIR, |
|
force_download: bool = False, |
|
width: int = 3840, |
|
height: int = 2400, |
|
filename: str | None = "warming.png", |
|
scale: int = 1, |
|
background: str = "#222", |
|
plot_background: str = "#202020", |
|
text_color: str = "#c9c9c9", |
|
axis_color: str = "#000", |
|
grid_color: str = "#000", |
|
) -> xy.Chart: |
|
rows = get_rows(url=url, force_download=force_download, data=data_dir) |
|
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=x_label, grid=show_x_grid), |
|
xy.y_axis( |
|
label=y_label, |
|
domain=(-6, 6), |
|
label_position="center", |
|
grid=show_y_grid, |
|
line=False, |
|
tick_values=[-5, -4, -2, 0, 2, 4, 5], |
|
), |
|
xy.theme( |
|
background=background, |
|
plot_background=plot_background, |
|
text_color=text_color, |
|
axis_color=axis_color, |
|
grid_color=grid_color, |
|
), |
|
title=title, |
|
width=width, |
|
height=height, |
|
) |
|
if filename: |
|
chart.to_png( |
|
filename, scale=scale |
|
) # native renderer: no browser / fonts needed |
|
print(f"Image {filename!r} saved.") |
|
return chart |
|
|
|
|
|
if __name__ == "__main__": |
|
import typer |
|
|
|
typer.run(get_chart) |
