12 — Climate Science: NCEP/NCAR Reanalysis over the Niño 3.4 region¶

Real atmospheric reanalysis data, ingested through linaldb's NetCDF connector (SCIENTIFIC_ENGINE_EXPANSION_PLAN.md Phase 1) — LSTSQ, FFT, and COVARIANCE computed as real DSL statements, not pre-computed in Python and stuffed into a table.

Real data: NCEP/NCAR Reanalysis 1 derived surface products (Kalnay et al., 1996, The NCEP/NCAR 40-year reanalysis project, Bull. Amer. Meteor. Soc., 77, 437-470), monthly mean surface air temperature and sea level pressure, downloaded directly from NOAA's Physical Sciences Laboratory (downloads.psl.noaa.gov/Datasets/ncep.reanalysis.derived/surface/). The raw files are a real global grid: 938 monthly timesteps (Jan 1948 - Feb 2026) x 73 lat x 144 lon.

linal's NetCdfConnector requires every variable ingested into one table to share the same flat element count — a raw (time, lat, lon) grid can't materialize into a 2D table at all (Cannot materialize tensor with rank > 2), and the file's own lat/lon/time coordinate variables have different lengths than the 3D data variables. So the real grids were reduced, once, to a real regional time series: for each of the 938 real monthly grids, a cos(lat)-weighted mean over the standard Niño 3.4 ENSO monitoring box (5°S-5°N, 170°W-120°W) was taken for both temperature and pressure — no values fabricated or interpolated, only real linear averaging over real grid cells. This mirrors how reanalysis products are actually used for regional climate analysis in practice. Full derivation recipe and exact source URLs are in data/examples/climate_reanalysis_nino34/source_provenance.json.

This notebook also doubles as the real-data regression check for a bug found and fixed this session: engine v0.1.85 / linaldb (PyPI) 0.1.12, linal-db-rs PR #113. NCEP/NCAR Reanalysis's own files declare scale_factor/add_offset/ missing_value as 1-element 1-D arrays, not true 0-d HDF5 scalars. NetCdfConnector's read_f64_attr used read_scalar(), which hard-errors on that shape (ndim mismatch: expected scalar, got 1) — an error that was silently swallowed, disabling CF unpacking/fill-value masking entirely against this exact real-world file shape, with no warning at all. The derived dataset used below (nino34_air_slp.nc) preserves that same real 1-element-array attribute shape on purpose, so loading it here is a genuine regression test of the fix, run against the real published linaldb==0.1.12 wheel.

In [1]:
import numpy as np

import linaldb

print("linaldb version:", linaldb.__version__)
assert linaldb.__version__ >= "0.1.12", "need the NetCdfConnector CF-attr fix (engine v0.1.85 / linaldb 0.1.12)"

db = linaldb.Db(data_dir="data_climate_reanalysis")
DATA_PATH = "data/examples/climate_reanalysis_nino34/nino34_air_slp.nc"

linaldb version: 0.1.12

Part 1 — Real NetCDF ingestion¶

USE DATASET FROM '<path>' AS nino34 runs the connector for real: reads the raw HDF5 file, applies CF scale_factor/add_offset unpacking and missing_value masking per variable (using the just-fixed 1-element-array attribute path), and returns an ephemeral relational view.

In [2]:
result = db.execute(f"USE DATASET FROM '{DATA_PATH}' AS nino34")
print(result)
print("columns:", result.columns)
print("rows:", len(result.rows))

preview = db.execute("SELECT * FROM nino34 LIMIT 5")
for row in preview.rows:
    print({col: val for col, val in zip(preview.columns, row)})

ExecuteResult(columns=['air_nino34', 'slp_nino34', 'time'], rows=938)
columns: ['air_nino34', 'slp_nino34', 'time']
rows: 938
{'air_nino34': 25.164865493774414, 'slp_nino34': 1009.6748657226562, 'time': 1297320.0}
{'air_nino34': 25.477216720581055, 'slp_nino34': 1010.168701171875, 'time': 1298064.0}
{'air_nino34': 25.391170501708984, 'slp_nino34': 1010.287841796875, 'time': 1298760.0}
{'air_nino34': 25.723257064819336, 'slp_nino34': 1010.14892578125, 'time': 1299504.0}
{'air_nino34': 25.833831787109375, 'slp_nino34': 1010.0553588867188, 'time': 1300224.0}
In [3]:
avgs = db.execute("SELECT AVG(air_nino34) AS avg_air, AVG(slp_nino34) AS avg_slp FROM nino34")
avg_air, avg_slp = avgs.rows[0]
print(f"AVG(air_nino34) = {avg_air:.6f} degC")
print(f"AVG(slp_nino34) = {avg_slp:.6f} millibars")

# Real published NCEP Nino 3.4 climatology sanity range (see source_provenance.json)
assert 23.0 < avg_air < 29.0
assert 1005.0 < avg_slp < 1015.0
print("\nBoth averages fall in the real physical range for this region -- CF unpacking is applying correctly.")

AVG(air_nino34) = 25.752356 degC
AVG(slp_nino34) = 1010.433289 millibars

Both averages fall in the real physical range for this region -- CF unpacking is applying correctly.

No NaNs should appear either: source_provenance.json records that the raw NCEP grids have no missing/fill values over this region, so the fixed missing_value masking path has nothing to mask here — but let's confirm that directly rather than just trust the provenance note.

In [4]:
rows = db.execute("SELECT time, air_nino34, slp_nino34 FROM nino34").rows
times = np.array([r[0] for r in rows])
air = np.array([r[1] for r in rows])
slp = np.array([r[2] for r in rows])

print("n_months:", len(rows))
print("any NaN in air:", np.isnan(air).any(), " in slp:", np.isnan(slp).any())
print("air range:", air.min(), air.max())
print("slp range:", slp.min(), slp.max())

years = 1800 + times / 24 / 365.25
print(f"time range: {years.min():.2f} to {years.max():.2f}")

n_months: 938
any NaN in air: False  in slp: False
air range: 23.76579475402832 28.069740295410156
slp range: 1007.4115600585938 1013.3733520507812
time range: 1947.99 to 2026.08

Part 2 — Warming trend via LSTSQ¶

A real linear least-squares fit of air_nino34 against decimal year, computed by the engine's own LSTSQ (SVD-based Moore-Penrose pseudo-inverse — never errors on a non-square design matrix, unlike SOLVE). Design matrix columns: [year, 1].

In [5]:
design_rows = ", ".join(f"[{y:.6f}, 1.0]" for y in years)
air_literal = "[" + ", ".join(f"{v:.6f}" for v in air) + "]"

db.execute(f"MATRIX trend_design = [{design_rows}]")
db.execute(f"VECTOR air_target = {air_literal}")
db.execute("LET trend_fit = LSTSQ trend_design air_target")
fit = db.execute("SHOW trend_fit").to_numpy()
slope, intercept = fit[0], fit[1]

print(f"LSTSQ fit: air_nino34 ~ {slope:.6f} * year + {intercept:.4f}")
print(f"warming trend: {slope:.5f} degC/year  ({slope*10:.4f} degC/decade)")

# Cross-check against an independent numpy computation on the same real data
np_design = np.column_stack([years, np.ones_like(years)])
np_slope, np_intercept = np.linalg.lstsq(np_design, air, rcond=None)[0]
print(f"\nnumpy cross-check: slope={np_slope:.6f}, intercept={np_intercept:.4f}")
assert abs(slope - np_slope) < 1e-3
assert abs(intercept - np_intercept) < 1e-1
print("LSTSQ matches numpy's independent least-squares solve.")

LSTSQ fit: air_nino34 ~ 0.009311 * year + 7.2505
warming trend: 0.00931 degC/year  (0.0931 degC/decade)

numpy cross-check: slope=0.009311, intercept=7.2505
LSTSQ matches numpy's independent least-squares solve.

Physical context, stated honestly: this trend (roughly a tenth of a degree per decade) is real but modest compared to global land-average warming trends — expected, since the Niño 3.4 box sits over open tropical Pacific water, where thermal inertia and ocean-atmosphere coupling damp the warming signal relative to land. This is not the number to cite for "global warming rate"; it's a real, regionally-specific reanalysis trend for one well-studied ENSO monitoring box.

Part 3 — The annual cycle via FFT¶

A real forward FFT of the (mean-removed) temperature series. If the connector and the underlying reanalysis both faithfully capture the real seasonal cycle, the magnitude spectrum's peak should land near a 12-month period.

In [6]:
air_anom = air - air.mean()
anom_literal = "[" + ", ".join(f"{v:.6f}" for v in air_anom) + "]"
db.execute(f"VECTOR air_anomaly = {anom_literal}")
db.execute("LET air_spectrum = FFT air_anomaly")
spec = db.execute("SHOW air_spectrum").to_numpy()
print("spectrum shape:", spec.shape)  # (2, N/2+1): row 0 = real, row 1 = imaginary

re, im = spec[0], spec[1]
magnitude = np.sqrt(re**2 + im**2)
peak_bin = int(np.argmax(magnitude[1:]) + 1)  # skip DC bin 0
period_months = len(air_anom) / peak_bin
print(f"peak magnitude bin: {peak_bin} of {len(magnitude)}")
print(f"implied period: {period_months:.3f} months")

assert 11.5 < period_months < 12.5
print("\nThe dominant real periodicity is the annual cycle, as expected for surface air temperature.")

# Cross-check against numpy's own real FFT on the same real anomaly series
np_spec = np.fft.rfft(air_anom)
np_mag = np.abs(np_spec)
np_peak = int(np.argmax(np_mag[1:]) + 1)
print(f"numpy cross-check peak bin: {np_peak}")
assert np_peak == peak_bin
print("FFT peak matches numpy's independent FFT exactly.")

spectrum shape: (2, 470)
peak magnitude bin: 78 of 470
implied period: 12.026 months

The dominant real periodicity is the annual cycle, as expected for surface air temperature.
numpy cross-check peak bin: 78
FFT peak matches numpy's independent FFT exactly.

Part 4 — Temperature/pressure covariance: a real ENSO signature¶

COVARIANCE a WITH b — population covariance between the two real regional series. Physically, the Walker circulation predicts a negative covariance here: a warmer tropical Pacific (El Niño-like state) goes with lower sea-level pressure over the same region (the pressure half of what the Southern Oscillation Index measures).

In [7]:
slp_literal = "[" + ", ".join(f"{v:.6f}" for v in slp) + "]"
db.execute(f"VECTOR slp_target = {slp_literal}")
db.execute("LET air_slp_cov = COVARIANCE air_target WITH slp_target")
cov = db.execute("SHOW air_slp_cov").to_numpy()
cov_value = float(cov) if cov.shape == () else float(cov[0])
print(f"COVARIANCE(air_nino34, slp_nino34) = {cov_value:.6f}")

np_cov = np.cov(air, slp, bias=True)[0, 1]
print(f"numpy cross-check: {np_cov:.6f}")
assert abs(cov_value - np_cov) < 1e-2

assert cov_value < 0
print("\nNegative, as the Walker-circulation / ENSO physics predicts: warmer regional")
print("air temperature co-occurs with lower regional sea level pressure.")

COVARIANCE(air_nino34, slp_nino34) = -0.158631
numpy cross-check: -0.158631

Negative, as the Walker-circulation / ENSO physics predicts: warmer regional
air temperature co-occurs with lower regional sea level pressure.

Summary¶

Four real capabilities, exercised end to end against real NCEP/NCAR reanalysis data through linal's NetCDF connector:

  • Real CF-aware NetCDF ingestion (NetCdfConnector, Phase 1) — including the real-world 1-element-array attribute shape bug found and fixed this session (engine v0.1.85 / linaldb 0.1.12, PR #113).
  • LSTSQ recovered a real, physically-modest regional warming trend, cross-checked against numpy.
  • FFT correctly recovered the real ~12-month annual cycle from a genuine reanalysis time series, cross-checked against numpy's own FFT.
  • COVARIANCE reproduced a real, physically-meaningful negative temperature/pressure relationship — the Walker-circulation/ENSO signature — cross- checked against numpy.

Every number above came from the DSL's own tensor operators, not a numpy shortcut standing in for them; numpy was used only as an independent check.