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.

In [1]:

runtime-identity

Copied!

# Record the exact code imported by this notebook kernel.
import cubedynamics as cd

print(cd.version_info())
# Record the exact code imported by this notebook kernel.
import cubedynamics as cd

print(cd.version_info())

CubeDynamics 0.1.0rc3
Artifact: development checkout
Code: git:a852137944cabf238e9321f2eecf09e483c9b9f5
Imported from: /home/runner/work/cubedynamics/cubedynamics/src/cubedynamics
Distribution: /opt/hostedtoolcache/Python/3.11.16/x64/lib/python3.11/site-packages

Define · Make the scientific contract visible in code¶

In [2]:

Copied!





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 [3]:

Copied!





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 [4]:

Copied!





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()

No description has been provided for this image

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. The first code cell prints cd.version_info() so a rendered result can be tied to a package path and, for development checkouts, a Git commit.