07 · Build a project-owned verb¶

Context¶

The core package supplies the grammar; projects supply domain verbs. A good custom verb states its input contract and returns labeled data.

Question¶

How much below-freezing exposure accumulated across the observed January minimum-temperature record?

Analysis story¶

We write one small verb, test direct and piped use, then visualize its state and magnitude outputs.

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.

Define · Make the scientific contract visible in code¶

In [1]:

import xarray as xr

def freezing_exposure(threshold=0.0):
    """Convert Celsius temperature to freezing state and magnitude."""
    def _op(cube):
        if cube.attrs.get("units") != "degC":
            raise ValueError("freezing_exposure requires units='degC'")
        if "time" not in cube.dims:
            raise ValueError("freezing_exposure requires a time dimension")

        state = cube <= threshold
        magnitude = (threshold - cube).where(state, 0.0)
        result = xr.Dataset({"state": state, "magnitude": magnitude})
        result.attrs.update(
            analysis="freezing_exposure",
            threshold_degC=float(threshold),
            source=cube.attrs.get("source", ""),
            is_synthetic=cube.attrs.get("is_synthetic", 0),
        )
        return result
    return _op

import xarray as xr def freezing_exposure(threshold=0.0): """Convert Celsius temperature to freezing state and magnitude.""" def _op(cube): if cube.attrs.get("units") != "degC": raise ValueError("freezing_exposure requires units='degC'") if "time" not in cube.dims: raise ValueError("freezing_exposure requires a time dimension") state = cube <= threshold magnitude = (threshold - cube).where(state, 0.0) result = xr.Dataset({"state": state, "magnitude": magnitude}) result.attrs.update( analysis="freezing_exposure", threshold_degC=float(threshold), source=cube.attrs.get("source", ""), is_synthetic=cube.attrs.get("is_synthetic", 0), ) return result return _op

Pipe · Test both forms on real observations¶

In [2]:

from pathlib import Path

import xarray as xr

# Find the repository from either a root-level documentation build or a kernel
# started beside this notebook, then open the checksum-controlled PRISM extract.
data_path = next(
    candidate / "tests" / "fixtures" / "real_data" / "prism_boulder_january_2024.nc"
    for candidate in (Path.cwd(), *Path.cwd().parents)
    if (candidate / "tests" / "fixtures" / "real_data" / "prism_boulder_january_2024.nc").exists()
)
prism = xr.open_dataset(data_path, engine="scipy").load()

# These assertions are part of the teaching contract: official source,
# canonical cube dimensions, complete daily time, and declared Celsius units.
assert prism.attrs["source"] == "PRISM Group, Oregon State University"
assert prism.attrs["is_synthetic"] == 0
assert prism.sizes == {"time": 30, "y": 24, "x": 24}
assert prism["tmax"].attrs["units"] == "degC"

import numpy as np
from cubedynamics import pipe

cube = prism["tmin"]
direct = freezing_exposure(-5.0)(cube)
through_grammar = (pipe(cube) | freezing_exposure(-5.0)).unwrap()

np.testing.assert_array_equal(through_grammar["state"], direct["state"])
np.testing.assert_allclose(through_grammar["magnitude"], direct["magnitude"])

from pathlib import Path import xarray as xr # Find the repository from either a root-level documentation build or a kernel # started beside this notebook, then open the checksum-controlled PRISM extract. data_path = next( candidate / "tests" / "fixtures" / "real_data" / "prism_boulder_january_2024.nc" for candidate in (Path.cwd(), *Path.cwd().parents) if (candidate / "tests" / "fixtures" / "real_data" / "prism_boulder_january_2024.nc").exists() ) prism = xr.open_dataset(data_path, engine="scipy").load() # These assertions are part of the teaching contract: official source, # canonical cube dimensions, complete daily time, and declared Celsius units. assert prism.attrs["source"] == "PRISM Group, Oregon State University" assert prism.attrs["is_synthetic"] == 0 assert prism.sizes == {"time": 30, "y": 24, "x": 24} assert prism["tmax"].attrs["units"] == "degC" import numpy as np from cubedynamics import pipe cube = prism["tmin"] direct = freezing_exposure(-5.0)(cube) through_grammar = (pipe(cube) | freezing_exposure(-5.0)).unwrap() np.testing.assert_array_equal(through_grammar["state"], direct["state"]) np.testing.assert_allclose(through_grammar["magnitude"], direct["magnitude"])

Figure · Interpret the custom verb's two outputs¶

In [3]:

import matplotlib.pyplot as plt

freezing_fraction = through_grammar["state"].mean(("y", "x"))
cumulative_exposure = through_grammar["magnitude"].sum("time")

fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.8), constrained_layout=True)
freezing_fraction.plot(ax=axes[0], marker="o", color="#365f79")
axes[0].set_title("Region below −5 °C each day")
axes[0].set_ylabel("Fraction of PRISM cells")
cumulative_exposure.plot(
    ax=axes[1], cmap="magma", cbar_kwargs={"label": "Degree-days below −5 °C"}
)
axes[1].set_title("Cumulative cold exposure")
plt.show()

import matplotlib.pyplot as plt freezing_fraction = through_grammar["state"].mean(("y", "x")) cumulative_exposure = through_grammar["magnitude"].sum("time") fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.8), constrained_layout=True) freezing_fraction.plot(ax=axes[0], marker="o", color="#365f79") axes[0].set_title("Region below −5 °C each day") axes[0].set_ylabel("Fraction of PRISM cells") cumulative_exposure.plot( ax=axes[1], cmap="magma", cbar_kwargs={"label": "Degree-days below −5 °C"} ) axes[1].set_title("Cumulative cold exposure") plt.show()

What the figure tells us¶

The custom verb separates occurrence from intensity. Its code remains a small project extension, while pipe supplies the shared composition language.

Try the next variation¶

Compose this verb with v.detect_events rather than embedding event detection.

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/custom_verb_project.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 ·

Learn the grammar · All vignettes