Companion notebook: Weather Data Lives on a Sphere¶
This notebook reproduces every figure and every benchmark number in the article Weather Data Lives on a Sphere. Our Statistics Should Too.
It is a thin, auditable layer over the repository's scripts: nothing here is reimplemented, each cell calls the same code the article used.
Quick mode (default). The two ERA5 cells reuse the cached netCDF files shipped
alongside this notebook (era5_heatdome_2021.nc, era5_jja_monthly_1991_2020.nc),
so no API key is needed. If the cached files are missing, the scripts fall back to
downloading from the Copernicus Climate Data Store, which requires a free CDS account
and an API token in ~/.cdsapirc (see the ERA5 section below).
Sections: 1. environment, 2. validation benchmarks, 3. property figures, 4. ERA5 applications.
1. Environment and versions¶
import sys
from importlib.metadata import version
print("python ", sys.version.split()[0])
for pkg in ["numpy", "geomstats", "matplotlib", "xarray", "cdsapi"]:
print(f"{pkg:<11}", version(pkg))
python 3.12.3 numpy 2.0.2 geomstats 2.8.0 matplotlib 3.11.0 xarray 2026.4.0 cdsapi 0.7.7
2. Validation benchmarks (the article's validation box)¶
Real weather has no planted ground truth, so correctness is checked on synthetic data with a known answer. Each line below states the property being tested and the measured result. These are the numbers quoted in the article.
2.1 Frechet mean (tracking engine)¶
B1 analytic ground truth, B2 agreement with GeomStats (after tightening its optimizer, whose default stops about 2.7 km short), B3 first-order optimality, B4 chordal-vs-geodesic gap growth, B5 date-line failure of the naive mean plus SO(3) equivariance, B6 planted trajectory recovery.
import tool1_centroid_tracker as tool1
tool1.bench()
B1 Recover a vMF blob's known center (analytic ground truth)
[PASS] frechet recovers center truth=(35.0, -100.0) got=(35.000,-100.000) err=0.00 km
B2 Hand-rolled Karcher == GeomStats FrechetMean
[PASS] karcher vs geomstats geodesic gap = 1.90e-04 km
B3 First-order optimality: ||sum w_i log_mu(p_i)|| ~ 0
[PASS] gradient residual ~ 0 residual = 5.99e-16
B4 Chordal (naive) vs geodesic (Frechet) gap grows with SPREAD
(a symmetric blob has gap 0 by symmetry - local flatness; the gap
appears only for asymmetric/spread data, as in figure 6 / hello_sphere)
tight cluster spread_radius= 3.1 deg chordal-vs-geodesic gap = 0.0 km
continental spread_radius= 33.4 deg chordal-vs-geodesic gap = 40.7 km
global spread spread_radius= 93.5 deg chordal-vs-geodesic gap = 7758.9 km
B5 Dateline + SO(3) equivariance
[PASS] frechet ok at dateline err = 0.00 km
[PASS] naive latlon FAILS err = 17805 km (avg lon -> 0.0 deg, wrong side)
[PASS] SO(3)-equivariant invariance error = 9.49e-05 km
B6 Recover a MOVING planted trajectory (10 days)
[PASS] trajectory recovered max_err=0.00 km total_path=5802 km
2.2 Principal Geodesic Analysis (variability engine)¶
B1 planted-axis recovery, B2 agreement with GeomStats TangentPCA, B3 no spurious axis
on isotropic data, B4 PGA versus naive lat/long PCA near a pole, B5 one-axis
reconstruction error, B6 SO(3) equivariance of the axes.
import tool3_geodesic_pca as tool3
tool3.bench()
B1 Recover a planted principal axis (ground truth) [PASS] axis recovered axis within 0.73 deg; captured 97.5% variance; mean off 69 km B2 Hand-rolled == GeomStats TangentPCA (integrity) [PASS] matches library axis angle diff 0.00e+00 deg; var-ratio diff 3.33e-16 B3 Isotropic -> no spurious axis; anisotropic -> real axis [PASS] isotropic ratio ~ 1 lambda1/lambda2 = 1.38 [PASS] anisotropic ratio high lambda1/lambda2 = 75.3 B4 PGA vs naive lat/long PCA near a pole (why geometry) [PASS] PGA axis correct at pole PGA off 0.8 deg [PASS] naive lat/long axis wrong naive off 39.2 deg (a 51x worse axis) B5 1-axis reconstruction error (justifies dimensionality reduction) [PASS] 1 axis reconstructs cloud RMS geodesic error 306 km (cloud spread 53 deg) B6 SO(3) equivariance (coordinate-free) [PASS] axes rotate with the data axis angle diff 0.00e+00 deg; var-ratio diff 4.44e-16
from IPython.display import Image, display
import spread_gap_figure
spread_gap_figure.make()
display(Image("spread_gap.png", width=700))
saved spread_gap.png tight cluster spread 3.1 deg gap 0.0 km continental spread 33.4 deg gap 40.7 km global spread spread 93.5 deg gap 7758.9 km
3.2 Algebra: the Frechet mean is SO(3)-equivariant, the naive mean is not (article Figure 2)¶
import equivariance_figure
equivariance_figure.make()
display(Image("equivariance.png", width=850))
left mean (lat,lon): Frechet (np.float64(22.0), np.float64(-100.6)) naive drift 15 km right mean (lat,lon): Frechet (86.0,-108.0) naive drift 782 km equivariance error |mu(R.X) - R.mu(X)| = 0.00013 km
saved equivariance.png
4. ERA5 applications¶
Data source. ERA5 is the ECMWF global reanalysis, distributed by the Copernicus Climate Data Store (CDS). The two cells below use:
- hourly 2 m temperature, June 24-30 2021, over a North America window
(
era5_heatdome_2021.nc, about 0.3 MB); - monthly-mean 2 m temperature, June-August 1991-2020, over a North America window
(
era5_jja_monthly_1991_2020.nc, about 4.5 MB).
If these cached files are present the cells run offline. Otherwise the scripts download
them once via cdsapi, which needs a free CDS account, a Personal Access Token in
~/.cdsapirc, and acceptance of the ERA5 licence on the dataset page. Credentials are
never stored in this repository.
4.1 Tracking the June 2021 heat dome (article Figure 3)¶
import era5_heatwave as heatwave
data = heatwave.load(heatwave.fetch()) # cached .nc if present
rows = heatwave.track(data)
heatwave.print_report(rows)
heatwave.plot(data, rows)
display(Image("era5_heatdome.png", width=650))
using cached era5_heatdome_2021.nc
day Frechet (lat,lon) naive (lat,lon) gap km 2021-06-24 ( 34.27,-103.57) ( 34.11,-103.51) 18.3 2021-06-25 ( 33.25,-106.78) ( 33.08,-106.82) 19.8 2021-06-26 ( 35.48,-112.04) ( 35.29,-112.35) 34.7 2021-06-27 ( 38.81,-116.67) ( 38.74,-116.97) 27.1 2021-06-28 ( 40.88,-117.52) ( 40.82,-117.80) 24.3 2021-06-29 ( 44.60,-117.72) ( 44.57,-117.88) 13.0 2021-06-30 ( 49.34,-114.99) ( 49.26,-114.67) 24.7 total Frechet centroid migration: 2635 km over 7 days
saved era5_heatdome.png
4.2 Interannual variability of summer heat, 1991-2020 (article Figure 4)¶
One anomaly heat-center per summer (that year's JJA mean minus the 30-year climatology), then PGA on the 30 centers.
import era5_pga_interannual as interannual
from tool3_geodesic_pca import GeodesicPCA
data30 = interannual.load(interannual.fetch()) # cached .nc if present
centers = interannual.yearly_centroids(data30)
pga = GeodesicPCA().fit(centers)
interannual.report(data30, centers, pga)
interannual.plot_pga(pga, centers, png="era5_pga_interannual.png",
title="Interannual variability of N. American summer heat (1991-2020)")
display(Image("era5_pga_interannual.png", width=650))
using cached era5_jja_monthly_1991_2020.nc
year heat-anomaly center (lat, lon) 1991 ( 52.67, -94.25) 1992 ( 34.22, -120.20) 1993 ( 38.47, -87.77) 1994 ( 40.17, -112.08) 1995 ( 45.90, -86.75) 1996 ( 44.77, -109.29) 1997 ( 36.31, -123.07) 1998 ( 48.67, -96.92) 1999 ( 52.05, -83.52) 2000 ( 41.45, -105.60) 2001 ( 48.89, -96.40) 2002 ( 43.14, -98.23) 2003 ( 48.15, -109.89) 2004 ( 48.90, -126.22) 2005 ( 46.80, -80.98) 2006 ( 45.99, -101.46) 2007 ( 44.65, -107.32) 2008 ( 56.16, -81.51) 2009 ( 41.69, -107.74) 2010 ( 41.14, -87.16) 2011 ( 34.60, -98.67) 2012 ( 44.31, -94.40) 2013 ( 49.73, -109.30) 2014 ( 47.50, -101.05) 2015 ( 41.92, -120.61) 2016 ( 40.68, -95.86) 2017 ( 44.57, -116.30) 2018 ( 39.68, -108.08) 2019 ( 41.22, -102.42) 2020 ( 41.10, -100.39) --- actionable result --- mean center of anomalous summer heat: (44.8, -102.5) axis 1 explains 77% of interannual variance (axis 2: 23%) axis 1 (+/-2 sigma) runs 47.1,-77.4 <-> 37.7,-123.9 (span ~3894 km)
saved era5_pga_interannual.png
That is the complete evidence base of the article: two benchmark suites, two property figures, two real-data applications. Methods are Karcher (1977) and Fletcher et al. (2004) as surveyed in Papillon, Sanborn, Mathe et al., Beyond Euclid (2025), computed with GeomStats.