04 · Read the analysis from left to right¶
Context¶
CubeDynamics makes the order of operations visible: a cube flows through small verbs instead of disappearing inside a long function call.
Question¶
Does the pipe grammar change the calculation, or only make the method easier to read and extend?
Analysis story¶
Follow observed temperature through selection, a local anomaly, a regional summary, and explicit output. Compare direct and piped standardization exactly. Each analytical cell displays its own result before the story continues.
Data used in this lesson¶
Every value comes from the PRISM Group at Oregon State University's AN91d daily 4 km climate product. This repository carries a small Boulder-region extract for 1–30 January 2024 so the lesson runs offline without replacing observations with generated values. The data validation page records source URLs, terms, checksums, bounds, units, and acceptance tests.
Prepare · Load checked real observations¶
REAL DATA · PRISM Boulder, January 2024. The shared setup verifies the fixture checksum before reading values.
from pathlib import Path
import hashlib
import json
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from IPython.display import display
from cubedynamics import pipe, verbs as v
# Run in the cloned repository or beside a downloaded notebook in the repo.
repo = next(p for p in (Path.cwd(), *Path.cwd().parents)
if (p / "tests/fixtures/real_data").is_dir())
def observed_cube(stem, variable):
path = repo / "tests/fixtures/real_data" / (stem + ".nc")
record = json.loads(path.with_suffix(".provenance.json").read_text())
assert hashlib.sha256(path.read_bytes()).hexdigest() == record["fixture_sha256"]
with xr.open_dataset(path, engine="scipy") as dataset:
assert not dataset.attrs["is_synthetic"]
result = dataset[variable].load() # Only this small, reviewed local extract.
result.attrs = {**dataset.attrs, **result.attrs}
assert result.dims == ("time", "y", "x")
assert np.all(np.diff(result.x) > 0) and np.all(np.diff(result.y) < 0)
assert bool(np.isfinite(result).all())
return result
cube = observed_cube("prism_boulder_january_2024", "tmax").rename("temperature")
assert cube.attrs["units"] == "degC"
plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 12, 'axes.titlesize': 13, 'axes.labelsize': 12, 'figure.facecolor': 'white', 'axes.facecolor': 'white', 'savefig.facecolor': 'white', 'figure.dpi': 140})
Pipe · Start with a temperature field¶
Where was the cold outbreak visible on 16 January?
field = cube.sel(time="2024-01-16")
fig, ax = plt.subplots(figsize=(5.4, 3.8), layout="constrained")
field.plot(ax=ax, cmap="magma", cbar_kwargs={"label": "Daily maximum (°C)"})
ax.set(title="PRISM · Boulder · 16 January 2024",
xlabel="Longitude (°E)", ylabel="Latitude (°N)")
plt.show()
Real PRISM daily maximum temperature, Boulder region, 16 January 2024. Selecting a date exposes one spatial face of the cube; north is up.
What changed?
Each cell is a gridded estimate, not a station reading. This map shows absolute temperature; it does not yet say how unusual the day was.
Pipe · Choose the time window¶
Which observations will define the short-period baseline?
window = (pipe(cube)
| v.apply(lambda c: c.sel(time=slice("2024-01-10", "2024-01-20")))).unwrap()
assert window.sizes["time"] == 11
fig, ax = plt.subplots(figsize=(5.4, 3.4), layout="constrained")
cube.isel(y=12, x=12).plot(ax=ax, color="0.65", label="Full fixture")
window.isel(y=12, x=12).plot(ax=ax, marker="o", color="#246b70", label="Selected days")
ax.set(title="One grid cell · select 10–20 January", xlabel="Date", ylabel="Daily maximum (°C)")
ax.legend(frameon=False)
plt.show()
PRISM maximum temperature at one Boulder grid cell. The apply/sel stage retains 10–20 January 2024 (colored markers) without changing their values.
What changed?
The grey observations remain in cube but not window. Subsequent anomalies use these eleven selected days, not the entire month or a climatology.
Pipe · Remove the local baseline¶
Where was 16 January colder than each pixel’s own selected-period mean?
departures = (pipe(window) | v.anomaly(dim="time")).unwrap()
fig, axes = plt.subplots(2, 1, figsize=(5.4, 6.5), layout="constrained")
window.sel(time="2024-01-16").plot(ax=axes[0], cmap="magma", cbar_kwargs={"label": "°C"})
departures.sel(time="2024-01-16").plot(ax=axes[1], cmap="RdBu_r", center=0,
cbar_kwargs={"label": "Departure (°C)"})
for ax, title in zip(axes, ["Before · absolute temperature", "After anomaly() · local departure"]):
ax.set(title=title, xlabel="Longitude (°E)", ylabel="Latitude (°N)")
plt.show()
PRISM, Boulder, 16 January 2024: absolute maximum temperature above and anomaly below. anomaly() subtracts each pixel’s 10–20 January mean; the diverging scale is centered on zero.
What changed?
Negative departures mean colder than that pixel’s baseline. The coldest absolute location need not have the largest departure. A short event-window baseline is not a long-term climate normal.
Pipe · Summarize region¶
When was the selected region cold relative to its local baselines?
regional_anomaly = (pipe(departures)
| v.mean(dim=("y", "x"), keep_dim=False)).unwrap()
assert regional_anomaly.dims == ("time",)
fig, ax = plt.subplots(figsize=(5.4, 3.4), layout="constrained")
regional_anomaly.plot(ax=ax, marker="o", color="#246b70")
ax.axhline(0, color="0.5", linewidth=0.8)
ax.set(title="After mean() · grid-cell average anomaly", xlabel="Date", ylabel="Departure (°C)")
plt.show()
The mean verb reduces PRISM y and x to a daily series for 10–20 January. The plotted quantity is the unweighted mean of grid-cell anomalies, in °C.
What changed?
Space has disappeared from the result, but time remains. This is an equal-cell average on a latitude/longitude grid, not an area-weighted regional estimate.
Pipe · Compare scales, not units¶
Does the pipe compute the same z-score as the explicit formula?
standardized = (pipe(window) | v.zscore(dim="time")).unwrap()
direct = (window - window.mean("time")) / window.std("time")
np.testing.assert_allclose(standardized, direct, rtol=1e-6, atol=1e-6)
# Standardization is dimensionless; replace the inherited temperature label.
standardized.attrs["units"] = "1"
fig, axes = plt.subplots(2, 1, figsize=(5.4, 5.8), layout="constrained")
axes[0].hist(window.values.ravel(), bins=24, color="#246b70")
axes[1].hist(standardized.values.ravel(), bins=24, color="#246b70")
axes[0].set(title="Before · selected PRISM values", xlabel="Daily maximum (°C)", ylabel="Cell-days")
axes[1].set(title="After zscore() · per-pixel scaling", xlabel="Standard deviations (unitless)", ylabel="Cell-days")
plt.show()
All PRISM grid-cell days in the selected window, before and after per-pixel zscore(). Histograms show distributions on different, explicitly labeled scales; the code asserts equivalence to the direct formula.
What changed?
Each pixel is centered and scaled by its own temporal variability. The pooled histogram is not fitted to a normal distribution, and its cell-days are not independent samples.
Pipe · Check the saved result¶
Does explicit NetCDF output preserve the analysis?
from tempfile import TemporaryDirectory
with TemporaryDirectory() as directory:
target = Path(directory) / "regional_anomaly.nc"
saved = (pipe(regional_anomaly) | v.to_netcdf(str(target), engine="scipy")).unwrap()
with xr.open_dataarray(target, engine="scipy") as reopened:
restored = reopened.load()
xr.testing.assert_identical(saved, restored)
output_table = pd.DataFrame({
"Check": ["Dimensions", "Observations", "Values and metadata"],
"Before": [str(saved.dims), str(saved.size), "Reference result"],
"After reopening": [str(restored.dims), str(restored.size), "Identical (asserted)"],
})
display(output_table)
| Check | Before | After reopening | |
|---|---|---|---|
| 0 | Dimensions | ('time',) | ('time',) |
| 1 | Observations | 11 | 11 |
| 2 | Values and metadata | Reference result | Identical (asserted) |
The PRISM regional-anomaly result is written with to_netcdf() and reopened. An identity assertion checks values, coordinates, names and attributes; the temporary example file is then removed.
What changed?
Export is an explicit side effect. The compact table is more useful here than another identical curve: it reports the tested round trip.
Figure sequence · Keep the method visible¶
Each result above comes from the code directly preceding it. The shared
definitions in scripts/visual_examples.py also generate Learn and reference
examples; the notebook is executed by the existing vignette runner.
What the figure tells us¶
The z-score assertion proves that the pipe is a composition language, not a new statistical definition. Selection defines the baseline, anomaly preserves space and time, and the spatial mean removes space. The export table verifies that saving the final series did not change its values or metadata.
Try the next variation¶
Change the selected dates and explain why the anomaly map changes even when the displayed day's absolute temperature does not.
Data used¶
| Field | Frozen analysis input |
|---|---|
| Provider | PRISM Group, Oregon State University |
| Product | AN91d daily 4 km time series |
| Dates | 2024-01-01 to 2024-01-30 |
| Fixture | tests/fixtures/real_data/prism_boulder_january_2024.nc |
| Provenance record | tests/fixtures/real_data/prism_boulder_january_2024.provenance.json |
The PRISM source reference describes current catalog support; the fixture record above identifies the observations used here. Data validation documents checksums and acceptance checks. The analytical baseline and thresholds belong to this story, not the provider.
Reproduce¶
Clone the repository, then run these commands from its root:
python -m pip install -e ".[vignettes]"
python scripts/run_vignettes.py docs/vignettes/grammar_basics.ipynb
No network is needed after installation. Open the downloaded notebook in Jupyter and run all cells to see the same figures. The website executes these cells during its strict build. Environment setup and the vignette contract explain the workflow.
See also¶
temperature · precipitation · anomaly · apply · mean · to_netcdf · zscore