Use cases across the stack

The library shape, as it ships¶

Figure 1 is the stack the cases run on, and Report 16 sets out its layering and ownership rules. One consequence of those rules shapes every case below. geocatalog and geopatcher are used on their own, next to the domain libraries, and the deployment infrastructure (YAML, hashing, the registry, artifacts) sits in pipekit-experiment rather than in pipekit core.

The running example: one methane plume¶

Every case processes the same kind of data, so their code lands on the same arrays.

The scene. An imaging spectrometer such as EMIT measures a radiance cube $\mat{Y} \in \mathbb{R}^{B \times H \times W}$, which has $B$ spectral bands over an $H \times W$ pixel grid with a pixel area of $\Delta A$ (60 m pixels give $\Delta A = 3600$ m²). We index the pixels by $n = 1, \ldots, HW$, and $\vect{y}_n \in \mathbb{R}^{B}$ is the spectrum of pixel $n$. A methane plume absorbs in the shortwave infrared, so it shows up as a small, spectrally specific dip in $\vect{y}_n$.

The enhancement. A matched filter estimates the column enhancement $\hat{\alpha}_n \in \mathbb{R}$ at each pixel from the background mean $\gvect{\mu}_y \in \mathbb{R}^{B}$ and covariance $\gmat{\Sigma}_y \in \mathbb{R}^{B \times B}$ of the spectra and the methane target signature $\vect{J}_h \in \mathbb{R}^{B}$. Its derivation and its limits are in radtrax · 1, stage 5.

The flux. Converting $\hat{\alpha}_n$ to a mass column $\Omega_n \in \mathbb{R}$ and masking the plume pixels $\mathcal{P} \subset \{1, \ldots, HW\}$ gives the integrated mass enhancement $\mathrm{IME} = \sum_{n \in \mathcal{P}} \Omega_n \Delta A$, and then the source rate $Q = U_\text{eff} \mathrm{IME} / \ell$ for an effective wind speed $U_\text{eff}$ and a plume length $\ell$ Varon et al., 2018. The derivation, and which length goes with which wind calibration, are in plumax · 1. Which package owns each step once radtrax and plumax exist is set by the GeoModels package boundaries. On this page, $\Omega$ is the methane mass column and not a spatial domain, and $\mathrm{IME}$, $Q$, $U_\text{eff}$ and $\ell$ are all scalars.

Every piece of this ships in geotoolz, and pipekit’s Graph composes them, because IMEEstimate takes two inputs (the mass column and the mask):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Shapes: B = 60 bands; H = W = 64 pixels; ΔA = 3600 m²
import pipekit_ops as pk
import geotoolz as gz

radiance = pk.Input("radiance")                          # Y: (B, H, W)          today
# α_n = J_hᵀ Σ_y⁻¹ (y_n − μ_y) / (J_hᵀ Σ_y⁻¹ J_h)
enh = gz.matched_filter.MatchedFilter(
    target=j_h, fit_on_call=True)(radiance)              # (B, H, W) → (H, W)   today
# Ω_n = α_n · ρ_CH₄    (ppm·m → kg m⁻²)
mass = gz.plume.ColumnToMass(gas="CH4", units_in="ppm_m",
                             units_out="kg_m2")(enh)     # (H, W) → (H, W)      today
# 𝒫 = {n : α_n > τ_Otsu}, connected components ≥ 20 px
mask = gz.plume.PlumeMask(threshold="otsu", min_area_px=20)(enh)  # bool (H, W) today
# IME = Σ_{n∈𝒫} Ω_n ΔA;   Q = U_eff · IME / ℓ
flux = gz.plume.IMEEstimate(wind_speed=3.0)(mass, mask)  # (H, W)² → dict       today

plume = pk.Graph(inputs={"radiance": radiance},
                 outputs={"enh": enh, "mask": mask, "flux": flux})   # today
out = plume(radiance=scene)
out["enh"].values.max()       # → 1518.8 ppm·m        (1500 injected)
out["flux"]
# → {'ime_kg': 595.1, 'length_m': 1994.5, 'wind_speed_m_s': 3.0,
#    'emission_rate_kg_s': 0.895, 'emission_rate_uncertainty_kg_s': 0.448}

plume (radiance in; enhancement, mask and flux out) is the object that most cases below pass around. Cases that need a single-input, single-output pipeline use a Sequential of the same operators instead.

What the cases need that does not ship yet¶

Running the original gallery against the pinned code leaves nine gaps (Table 1). Each gap gets a letter, and every proposed (X) marker in the cases points at one of them:

Exploration and batch processing¶

Case 1 · Notebook exploration (pipekit → geotoolz)¶

TL;DR. Look at a scene through a cleaned, clipped index, and print what each step does to it.

Problem. A Sentinel-2 scene $\mat{X} \in \mathbb{R}^{B \times H \times W}$ with a QA60 band. We want to mask the clouds, stretch each band, compute the NDVI and look at the map.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Shapes: B bands (QA60 last), H × W pixels
import numpy as np
import pipekit_ops as pk
import geotoolz as gz

viz = pk.Sequential([
    # x̃ = NaN where QA60 bits 10, 11 are set
    gz.mask.ApplyMask(mask=gz.qa.MaskClouds(qa_band="QA60", bits=[10, 11])),
                                                    # (B, H, W) → (B, H, W)  today
    pk.Tap(lambda gt: print(f"NaN: {np.isnan(gt.values).mean():.1%}")),  # today
    # x̂_b = (x̃_b − p₂)/(p₉₈ − p₂), clipped to [0, 1]
    gz.radiometry.PercentileClip(lower=2.0, upper=98.0),  # (B, H, W)   today
    # NDVI = (ρ_nir − ρ_red)/(ρ_nir + ρ_red + tol)
    gz.indices.NDVI(nir=7, red=3),                  # (B, H, W) → (1, H, W)  today
])
ndvi = viz(scene)                                   # GeoTensor          today

Library map. pipekit and geotoolz, and the scene came from an earlier georeader call. Notes. The case is unchanged from the original gallery apart from the names: masks live in geotoolz.qa and are applied with geotoolz.mask.ApplyMask, and PercentileClip takes lower and upper.

Case 2 · ETL pipelines (geocatalog → pipekit → geotoolz → georeader)¶

TL;DR. Run one saved pipeline over every scene in a catalog, in parallel, and write a COG per scene.

Problem. A catalog of $N_\text{scene}$ scenes and a pipeline $f$ defined in YAML. We compute $f(\mat{X}_m)$ for $m = 1, \ldots, N_\text{scene}$ and write each result, while tolerating individual failures.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Shapes: N_scene scenes; each X_m: (B, H, W)
from geocatalog import open_catalog
from geocatalog.load import load_raster
from pipekit_ops import ProcessMap
from pipekit_experiment.adapters.hydra import HydraConfigLoader

catalog = open_catalog("s3://bucket/s2_eu.parquet")      # N_scene rows, lazy today
etl = HydraConfigLoader.from_yaml(open("ndvi_etl.yaml").read())  # nested     today

# f(X_m) for m = 1 … N_scene, on 8 processes; failures logged, not raised
CatalogPipeline(                                         # proposed (C)
    catalog=catalog,
    op=etl,
    loader=load_raster,                          # (catalog, slice) → GeoTensor  today
    parallel=ProcessMap(etl, n_workers=8, on_error="log_and_continue"),  # today
).run()                                                  # arguments sketched

Library map. geocatalog (iteration and loading), pipekit (parallelism), pipekit-experiment (the YAML loader) and geotoolz (the operators named in the YAML, ending in geotoolz.io.WriteCOG), with georeader underneath the loader. Notes. The YAML is the source of truth. A hand-written nested YAML with _target_ keys loads today, but dumping a built nested pipeline does not round-trip (proposal A). Every step must pickle for ProcessMap, and pipekit_ops.check_pickleable(etl) lists the operators that will not (Report 2).

Case 8 · Data validation as operators (pipekit → pipekit-array → geotoolz)¶

TL;DR. Put permanent runtime checks into the pipeline, so that CI and production both fail loudly on bad data.

Problem. Before and after a correction step, we check the carrier’s shape, CRS, resolution, band schema, valid fraction and value range.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# Shapes: X: (B, H, W)
import pipekit_ops as pk
import geotoolz as gz

qc = pk.Sequential([
    pk.AssertShape((5, None, None)),                    # B = 5, any H, W    today
    pk.AssertHasAttribute("crs"),                       # today
    AssertCRSEquals("EPSG:32630"),                      # proposed (E)
    AssertResolutionWithin((9.5, 10.5)),                # proposed (E)
    # valid fraction = |{n : X_n finite}| / (H W) ≥ 0.5
    AssertValidFraction(min_valid=0.5),                 # proposed (D)
    AssertValueRange(min_val=0, max_val=10_000, on_fail="warn"),  # proposed (D)
])
etl = pk.Sequential([
    gz.mask.ApplyMask(mask=gz.qa.MaskClouds(qa_band="QA60", bits=[10, 11])),  # today
    qc,                                                 # raises QCError     today
    gz.radiometry.SimpleAtmosphericCorrection(),        # today; arguments sketched
    qc,
    gz.indices.NDVI(nir=3, red=2),                      # today
])

Notes. The split into three tiers (Table 2) is still the right answer to “where does a check go?”, but only the first tier exists. Until D lands, pk.AssertCallable(lambda gt: np.isfinite(gt.values).mean() >= 0.5) covers the valid fraction, at the cost of a closure (forbid_in_yaml). There is no warn mode in pipekit’s asserts, which raise QCError, so Quarantine is the non-raising alternative. TOAToBOA from the original gallery does not exist, and SimpleAtmosphericCorrection is the nearest operator to it.

Case 11 · Cross-sensor fusion (the Graph case)¶

TL;DR. Combine optical, SAR and a DEM on one grid and classify the stack, keeping the fused stack as a second output.

Problem. The inputs are $\mat{X}_\text{opt} \in \mathbb{R}^{B_o \times H \times W}$, $\mat{X}_\text{sar} \in \mathbb{R}^{B_s \times H' \times W'}$ and $\mat{X}_\text{dem} \in \mathbb{R}^{H'' \times W''}$ on three grids. We align them to the optical grid, stack them and predict a class per pixel.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# Shapes: B_o optical, B_s SAR bands; H × W the optical grid; C classes
import pipekit_ops as pk
import geotoolz as gz

opt, sar, dem = pk.Input("optical"), pk.Input("sar"), pk.Input("dem")   # today
opt_c = gz.mask.ApplyMask(mask=gz.qa.MaskClouds(qa_band="QA60", bits=[10]))(opt)
                                                    # (B_o, H, W)             today
# σ⁰_dB = 10 log₁₀ σ⁰, after a 7 × 7 Lee filter
sar_db = LinearToDB()(gz.restore.DespeckleLee(window=7)(sar))
                                         # (B_s, H', W')  today + proposed (I)
fused = AlignAndStack(target="optical")(opt_c, sar_db, dem)
                                         # → (B_o + B_s + 1, H, W)  proposed (I)
crops = gz.learn.ModelOp(model=classifier)(fused)   # → (C, H, W)             today

g = pk.Graph(inputs={"optical": opt, "sar": sar, "dem": dem},
             outputs={"crop_map": crops, "fused": fused})    # today

Notes. The Graph topology is carrier-agnostic and works today, and the running example is a smaller instance of the same idea. Choosing the target grid decides what information is resampled away, so it stays a domain decision in geotoolz. The original gallery said that Graph.summary() catches a missing input at construction. However, it does not exist, and Graph checks for missing inputs when it is called (Report 2).

Machine learning¶

Case 3 · Training and sliding-window inference (geocatalog → pipekit-train → geopatcher)¶

TL;DR. Train on random chips from a catalog, then run the model over a whole scene in overlapping patches, with one preprocessing object on both paths.

Problem. Training draws chips $\mat{X}_k \in \mathbb{R}^{B \times 256 \times 256}$ for $k = 1, \ldots, K$. Inference splits a scene into $N_\text{patch}$ patches at a stride of 224 (32 pixels of overlap), predicts each patch $p$, and stitches the predictions by overlap-add (Report 7).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
# Shapes: B bands; chips 256 × 256; K training chips; N_patch inference patches
import pipekit_ops as pk
import geotoolz as gz
import geopatcher as gp
from geopatcher import spatial
from geopatcher.integrations.pipekit import ApplyToChips, GridSampler, MergePatches
from pipekit_train import IterableDataset, TrainingLoop

preprocess = pk.Sequential([                         # shared by both paths
    gz.radiometry.ToFloat32(),                # (B, H_win, W_win) → (B, H_win, W_win)  today
    gz.radiometry.PercentileClip(lower=2.0, upper=98.0),        # today
])
augment = pk.Sequential([
    gz.augment.RandomFlip(seed=0),                   # today
    gz.augment.RandomRotate90(seed=0),               # today
])

# Training: (augment ∘ preprocess)(X_k), y_k for k = 1 … K
train = IterableDataset(
    source=lambda: ((augment(preprocess(x)), y) for x, y in chips(catalog)),
    content_hash=catalog_hash)                       # today; chips() is user code
model, _ = TrainingLoop(model_op=gz.learn.ModelOp(model=net), dataset=train,
                        max_steps=10_000, batch_size=32).run()  # today; args sketched

# Inference: ŷ_n = Σ_p w_n^(p) ŷ_n^(p) / Σ_p w_n^(p)
patcher = gp.SpatialPatcher(
    geometry=spatial.geometry.Rectangular(size=(256, 256)),     # today
    sampler=spatial.sampler.RegularStride(step=(224, 224)),     # today
    window=spatial.window.Hann(),                               # today
    aggregation=spatial.aggregation.OverlapAdd(),               # today
)
infer = pk.Sequential([
    GridSampler(patcher=patcher),                    # field → N_patch patches  today
    ApplyToChips(operator=preprocess | model),   # (B, 256²) → (1, 256²)  today
    MergePatches(aggregation=spatial.aggregation.OverlapAdd(),
                 domain=field.domain),               # N_patch patches → (1, H, W)  today
])
prediction = infer(gp.RasterField(scene))           # today

Library map. geocatalog (chip sources), pipekit (composition), pipekit-train (the loop and the dataset), geopatcher with its [pipekit] extra (patching), geotoolz (operators and ModelOp) and georeader (the field). Notes. Using the same preprocess object on both paths is the guarantee against train/serve skew, and it holds. Two names moved since the original gallery:

• pa.ModelOp became gz.learn.ModelOp, since pipekit-array has no ModelOp.

• Stitch became MergePatches.

The original also had catalog.iter_random_slices, which does not exist. Random chips come from geopatcher’s spatial.sampler.Random over a RasterField, or from your own iterator, as here.

Case 10 · Active learning (geopatcher → pipekit → geotoolz)¶

TL;DR. Score every chip by the model’s uncertainty and label the hundred most uncertain ones first.

Problem. For the patches $p = 1, \ldots, N_\text{patch}$, we compute a per-patch score $\bar{\sigma}^{(p)} \in \mathbb{R}_{\ge 0}$, which is the mean predictive standard deviation over the patch’s pixels $\mathcal{N}^{(p)}$, and we select the 100 patches with the largest $\bar{\sigma}^{(p)}$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# Shapes: N_patch patches of 256 × 256; 100 selected
import heapq
import numpy as np
import pipekit_ops as pk
import geopatcher as gp
from geopatcher import spatial

score = pk.Sequential([
    preprocess,                                      # today
    gz.learn.ModelOp(model=uq_model, method="predict_std"),  # → (1, 256, 256)  today
    # σ̄^(p) = (1 / |𝒩^(p)|) Σ_{n∈𝒩^(p)} σ_n^(p)
    pk.Lambda(lambda s: float(np.nanmean(s)), name="mean"),  # → scalar     today
])
patcher = gp.SpatialPatcher(
    geometry=spatial.geometry.Rectangular(size=(256, 256)),
    sampler=spatial.sampler.RegularStride(step=(256, 256)),   # no overlap
    window=spatial.window.Boxcar(),
    aggregation=spatial.aggregation.ByIndex(),   # per patch, no stitch  today
)
# the 100 largest σ̄^(p): the patches to label next
top = heapq.nlargest(100, patcher.split(field),
                     key=lambda p: score(p.data))    # user; or proposed (H)

Notes. ByIndex is right because the scores stay per patch rather than being stitched back into a field. The original used pa.MeanScalar(field="sigma"). However, pipekit-array’s MeanScalar reduces over array axes (axis=) and has no field argument, so on a plain array pipekit_array.MeanScalar() works too. Scoring is the expensive part, which makes geopatcher.run.parallel_map the place to parallelise it.

Serving¶

Case 4 · A LitServe inference API (georeader → pipekit → geotoolz)¶

TL;DR. Serve the plume pipeline behind a GPU endpoint that reads the requested scene and returns the flux with the pipeline’s hash.

Problem. A request carries a scene URL and a bounding box. We read that window, run the pipeline and return the result together with an identifier that pins exactly which pipeline produced it.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Shapes: one request → one window X: (B, H_win, W_win)
import litserve as ls
import tempfile
from georeader.rasterio_reader import RasterioReader
from georeader.read import read_from_bounds
from pipekit_experiment import LocalModelRegistry
from pipekit_experiment.adapters.hydra import HydraConfigLoader

class PlumeAPI(ls.LitAPI):                            # user
    def setup(self, device):
        self.pipe = HydraConfigLoader.from_yaml(open("plume_v3.yaml").read())  # today
        # hash = SHA-256(canonical config): stable across processes
        self.hash = LocalModelRegistry(tempfile.mkdtemp()).store(self.pipe)    # today

    def decode_request(self, req):
        return read_from_bounds(RasterioReader(req["url"]), tuple(req["bbox"]))
                                                     # → GeoTensor (B, H_win, W_win)  today

    def predict(self, scene):
        return self.pipe(scene)                      # → flux dict             today

    def encode_response(self, out):
        return {**out, "pipeline_hash": self.hash}   # user

ls.LitServer(PlumeAPI(), accelerator="gpu").run(port=8000)   # user

Notes. The original’s georeader.read_from_url and pk.config_hash do not exist. Reading is read_from_bounds on a RasterioReader (which opens URLs through GDAL), and the hash is the registry’s content hash. Storing the same pipeline twice gives the same 64-character SHA-256, and changing one parameter changes it (we checked this). Proposal B would expose the hash without the throwaway registry. Heavy state lives in setup and never in predict, which is operational discipline rather than framework.

Case 5 · A multi-pipeline backend (FastAPI)¶

TL;DR. Serve several named pipelines, each self-describing and pinned by hash, and let a promotion swap which one a name points at.

Problem. A map from names to pipelines. Each endpoint reads, runs and writes a COG, and the API lists every pipeline’s config and hash.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# Shapes: one request → X: (B, H_win, W_win) → COG
from fastapi import FastAPI
from georeader.save import save_cog
from pipekit_experiment import LocalModelRegistry

app = FastAPI()                                      # user
reg = LocalModelRegistry("registry/")                # today
# name ↦ hash, e.g. "plume" ↦ "9246b8…"; promotion = re-pointing the tag
for name in ("plume", "ndvi_etl"):
    PIPES[name] = reg.load(reg.resolve_tag(name))    # today (flat pipelines)

@app.post("/process/{name}")
def process(name: str, req: Request):                # user
    out = PIPES[name](read_window(req))              # read_window: case 4  today
    save_cog(out["enh"], f"s3://outputs/{req.id}.tif")   # → COG            today
    return {"cog": f"s3://outputs/{req.id}.tif", "hash": reg.resolve_tag(name)}

@app.get("/pipelines/{name}")
def describe(name: str):
    return PIPES[name].get_config()                  # self-documenting      today

Notes. ConfigMixin makes get_config() work for every operator, so the API documents itself. The registry’s tags (tag, resolve_tag) are the promotion mechanism: we point plume at a new hash, and every server that resolves the tag picks it up. There is one limit, which is that registry.load rebuilds flat pipelines only, so a nested pipeline is stored for its hash and loaded from its YAML (proposal A closes the gap).

Case 6 · User-uploaded pipelines, sandboxed¶

TL;DR. Let users paste a pipeline and run it, but only with operators on an allowlist.

Problem. Untrusted text in and a pipeline out, with every class it names checked against an allowlist before anything is instantiated.

1
2
3
4
5
6
7
8
9
10
11
12
13
import pipekit_ops as pk
import geotoolz as gz

# 𝒜 = {operator classes users may instantiate}
for cls in (gz.indices.NDVI, gz.radiometry.PercentileClip):
    pk.register(cls)                                 # today
pk.allowed()                                         # → ('geotoolz…NDVI', …)  today

try:
    pipe = pk.loads_sandboxed(user_text)             # flat JSON           today
    # pipe = loads_yaml_sandboxed(user_yaml)         # nested YAML      proposed (F)
except PermissionError as err:                       # not in 𝒜            today
    show_error(err)                                  # user

Notes. The allowlist is the right shape, since loads_sandboxed refuses any class that is not registered. There are two gaps:

• Flat only. The JSON form is flat, so a user pipeline is one operator until F walks a nested YAML and checks each _target_.

• No introspection. The original’s pk.registry.public_operators(...) and UntrustedTargetError do not exist. Instead, you register classes explicitly, which is arguably safer.

forbid_in_yaml operators (closures) can never be loaded this way, which is the point. Resource limits are still the deployment’s job.

Case 7 · A tile server (z/x/y → PNG)¶

TL;DR. Render map tiles of a derived layer on demand.

Problem. For the tile (z, x, y), we read the matching 256 × 256 window in Web Mercator, run the layer’s pipeline and return a PNG.

1
2
3
4
5
6
7
8
9
10
# Shapes: tile (3, 256, 256) uint8
import mercantile
from georeader.read import read_from_bounds

@app.get("/tiles/{layer}/{z}/{x}/{y}.png")           # user (FastAPI)
def tile(layer: str, z: int, x: int, y: int):
    bounds = mercantile.xy_bounds(x, y, z)           # user
    gt = read_from_bounds(reader, bounds, crs_bounds="EPSG:3857")  # (B, H_win, W_win) today
    rgb = LAYERS[layer](gt)               # e.g. TrueColor | StretchToUint8  today
    return png_response(rgb.values)                  # user

Notes. geotoolz’s viz.TrueColor and viz.StretchToUint8 make the layer pipeline, whereas the original’s to_png_bytes and read_tile helpers are not stack API. For an indexed archive, geocatalog finds the files under the tile and geocatalog.load.load_raster mosaics them. A tile cache in front matters more than Python speed.

Reproducibility and operations¶

Case 9 · A pinned, hashed regulatory artifact (pipekit-experiment)¶

TL;DR. Freeze exactly how a published number was produced, and be able to rerun it years later.

Problem. We record the pipeline, the hashes of the models it uses, the dependency lock and the input identifiers in one file that can be checked and replayed.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from pipekit_experiment import InferenceArtifact, LocalModelRegistry
from pipekit_experiment.adapters.hydra import HydraConfigLoader

# At publication
yaml_text = open("plume_v3.yaml").read()
art = InferenceArtifact(                             # today
    pipeline_yaml=yaml_text,
    pinned_model_hashes={"classifier": reg.store(model, weights=w)},  # today
    backend_info={"python": "3.12", "image": "sha256:…"},
    deps_lock=open("pixi.lock").read(),
    metadata={"scene": "EMIT_L1B_RAD_…", "scene_sha256": "ab12…", "doi": "10.xxxx/…"},
)
art.save("artifacts/plume_2024Q3.json")              # today

# Years later
art = InferenceArtifact.load("artifacts/plume_2024Q3.json")      # today
pipe = HydraConfigLoader.from_yaml(art.pipeline_yaml)            # today
models = art.reload_models(reg)                      # today; weights round-trip

Notes. The original gallery put this in pipekit core as pk.repro. It shipped in pipekit-experiment instead, along with TrainingArtifact for the training side. Bit-exact reruns still need a container image hash, which is why backend_info exists. The operators inside can come from any library, because each one records its class path.

Case 12 · Orchestrator task units (Prefect, Dagster)¶

TL;DR. Hand each scene to an orchestrator task that retries safely and writes to a path that depends only on its inputs.

Problem. One task per slice. A retry must produce the same output path, so the key is a function of (pipeline hash, slice) and is never random.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from prefect import flow, task                       # user
from geocatalog import open_catalog
from geocatalog.load import load_raster
from georeader.save import save_cog

@task(retries=3, retry_delay_seconds=60)
def run_one(pipe_yaml: str, sl, catalog_uri: str) -> str:
    pipe = HydraConfigLoader.from_yaml(pipe_yaml)    # today
    cat = open_catalog(catalog_uri)                  # today
    # key = hash(pipeline) ‖ hash(slice): identical on every retry
    key = f"{pipeline_hash(pipe)[:12]}-{hash(sl):x}"  # proposed (B)
    out = pipe(load_raster(cat, sl))                 # → GeoTensor             today
    save_cog(out, f"s3://daily/{key}.tif")           # idempotent by key       today
    return key

@flow
def daily(date: str):
    cat = open_catalog("s3://…/catalog.parquet")
    for sl in cat.iter_slices(resolution=(60.0, 60.0)):   # GeoSlice/row  today
        run_one.submit(PIPE_YAML, sl, "s3://…/catalog.parquet")

Notes. GeoSlice is a frozen dataclass, so it pickles to a worker and hashes stably. The original used GeoSlice.from_dict and iter_slices_for_date, neither of which exists, so we filter the catalog with query(catalog, slice_) for a date window first. The orchestrator owns retries and fan-out, whereas the stack owns what each task computes.

Case 13 · Pipeline diffing and A/B regression¶

TL;DR. Before promoting a new pipeline version, show what changed in its config and how much its outputs moved.

Problem. Two pipelines, $f$ (current) and $f'$ (candidate). We report the config difference and, for a sample of $N_\text{scene}$ scenes, the per-scene $\mathrm{RMSE}_m = \Vert f(\mat{X}_m) - f'(\mat{X}_m) \Vert_2 / \sqrt{HW}$ and its distribution.

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: N_scene sample scenes; outputs (H, W)
import numpy as np
from pipekit_experiment.adapters.hydra import HydraConfigLoader as H

old = H.from_yaml(open("plume_v3.yaml").read())      # today
new = H.from_yaml(open("plume_v4.yaml").read())      # today
diff = config_diff(old.get_config(), new.get_config())   # proposed (G)

# RMSE_m = ‖f(X_m) − f′(X_m)‖₂ / √(HW),  m = 1 … N_scene
def enh(pipe, x):
    return pipe(x)["enh"].values                    # (H, W)               today
rmse = np.array([np.sqrt(np.nanmean((enh(old, x) - enh(new, x)) ** 2))
                 for x in sample_scenes])            # (N_scene,)           user

Notes. A mean RMSE hides a bimodal regression, so we look at the whole $\{\mathrm{RMSE}_m\}$ distribution or, better, at the change in the quantity that matters ($\Delta Q_m$ for the plume). For xarray outputs, xrtoolz.metrics.RMSE(variable, dims) computes it over named dimensions. pipekit-evaluate’s Benchmark is the structured version of this question for estimators against references, and extending it to “pipeline A against pipeline B” is proposal G.

New cases from the geoapps design¶

The geoapps platform (validation, attribution, watchlists, research) runs on this stack, and its design plan adds six shapes that the original gallery did not have. They matter because they stress the stack where real platforms live: versioned reruns, human verdicts and models that are retrained from them. The app-side pieces (the database, the queue, the UI) belong to geoapps, so the cases below show only the stack’s half.

Case 14 · A validation queue with priorities (pipekit → geotoolz)¶

TL;DR. Rank new predicted plumes so that analysts validate the most urgent ones first, and nothing unvalidated goes downstream.

Problem. Each predicted detection $d = 1, \ldots, N_\text{det}$ has a flux estimate $\hat{Q}_d \in \mathbb{R}_{\ge 0}$, an age $a_d \in \mathbb{R}_{\ge 0}$ and a flag $\delta_d \in \{0, 1\}$ (would validating it trigger a notification?). The priority $\pi_d \in \mathbb{R}_{\ge 0}$ combines them,

with the age weight

where $g$ increases with flux, $a_\text{notify} \in \mathbb{R}_{>0}$ is the notification window and $a_\text{decay} \in \mathbb{R}_{>0}$ is the decay after it. The code below computes (1) from the running example’s evidence.

1
2
3
4
5
6
# Shapes: N_det detections; each with (H, W) enhancement and mask
evidence = plume(radiance=scene)                     # running example     today
q_hat = evidence["flux"]["emission_rate_kg_s"]       # Q̂_d (kg s⁻¹)        today
sigma = evidence["flux"]["emission_rate_uncertainty_kg_s"]   # σ_Q        today
# π_d = g(Q̂_d) (1 + δ_d) φ(a_d)
priority = g(q_hat) * (1 + would_notify) * phi(age)  # geoapps-core       user

Notes. The stack’s job ends at the evidence: an enhancement map, a mask, $\hat{Q}_d$ and its uncertainty $\sigma_Q$, all tied to the pipeline hash that made them. The platform’s rule that a human has the final say is enforced in geoapps’ database, not in the stack. Each verdict becomes a label anchored to (scene, geometry), which case 18 consumes.

Case 15 · Scheduled pipeline families: discovery, monitoring, reanalysis, backfill¶

TL;DR. Run the same steps on four schedules, with every output keyed so that reruns are no-ops and a crashed backfill resumes.

Problem. A pipeline $f$ with a version tuple $v$ (pipeline hash, model hashes, input product versions). Each step’s output is keyed by $\mathrm{key} = (\text{step}, \text{input id}, v)$. The families differ only in when they run, what they cover and where they write (Table 3):

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: N_scene new scenes per poll
from geocatalog.sources import STACSource
from geocatalog.storage import CatalogBundle

# discovery: new items in an AOI and window, recorded with their query
bundle.ingest(STACSource("https://planetarycomputer…/api/stac/v1"),
              bounds=aoi, interval=window, collection="emit-l1b-rad")   # today
# key = (step, input id, v): rerun ⇒ no-op
v = (registry.store(plume_pipe), model_hash, product_version)   # today (registry)
for sl in bundle.catalog.iter_slices(resolution=(60.0, 60.0)):  # today
    key = (step, sl, v)                              # GeoSlice hashes stably  today
    if store.has(key): continue                      # geoapps job queue     user
    store.put(key, plume_pipe(load_raster(bundle.catalog, sl)))   # today

Notes. The stack supplies three things the families need:

• Discovery with provenance. CatalogBundle.ingest records the query that found each item, so a reanalysis can replay discovery (Report 6).

• Stable keys. GeoSlice hashes, and the registry hashes pipelines and models.

• Resumable patch runs. geopatcher.observe.PatchJournal resumes a killed run inside one large scene.

Monitoring’s non-detections (“the site was observable and nothing was seen”) are what make a persistence estimate possible, and without them it is just a detection count. Reanalysis never overwrites, since it writes its own stream under one frozen $v$.

Case 16 · What-if reruns in a sandbox¶

TL;DR. Rerun a published detection with one assumption changed (such as the wind product), and compare the two side by side without touching production.

Problem. We swap $U_\text{eff}$ (ERA5) for $U'_\text{eff}$ (GEOS-FP) with the mask held fixed. Since $Q = U_\text{eff} \mathrm{IME} / \ell$ (plumax · 1), the flux scales exactly:

1
2
3
4
5
# Shapes: (H, W) enhancement and mask, reused
# Q′/Q = U′_eff / U_eff, with 𝒫 fixed
q = {u: gz.plume.IMEEstimate(wind_speed=u)(mass_gt, mask_gt)["emission_rate_kg_s"]
     for u in (3.0, 4.5)}                            # today
q[4.5] / q[3.0]                    # → 1.5000000000000002   (= 4.5 / 3.0)

Notes. That ratio is an exact test oracle (gate J4 in the geoapps plan), and it holds on the shipped IMEEstimate. The stack makes what-ifs cheap because a pipeline is data, so we change one constructor argument, rehash, and the result lands under a new version tuple in a sandbox/<user>/<name> stream. The attribution can still flip under the new wind direction, which is the interesting part for analysts.

Case 17 · Catalog explorer: preview, then ingest¶

TL;DR. Let an analyst browse a public catalog, preview a few scenes over their area, and ingest the collection only if it is worth it.

Problem. For an AOI and a time window, we list the matching items, load a cheap preview of one, and either stage the collection or discard it.

1
2
3
4
5
6
7
8
9
10
11
12
13
from geocatalog import GeoSlice, query
from geocatalog.load import load_raster
from geocatalog.sources import STACSource
from geocatalog.staging import stage

src = STACSource("https://earth-search.aws.element84.com/v1")   # today
rows = list(src.query(bounds=aoi, interval=window,
                      collection="sentinel-2-l2a", limit=20))   # 20 rows  today
# preview: one scene at 100 m instead of 10 m (100× fewer pixels)
preview = load_raster(catalog_of(rows), GeoSlice(bounds=aoi, interval=window,
                      resolution=(100.0, 100.0), crs="EPSG:32630"))   # today
if analyst_says_yes(preview):                        # geoapps UI           user
    stage(catalog_of(rows), dest="staging/", parallel=8)   # download + cache   today

Notes. “ETLs are user-callable” is the geoapps rule that drives this case. The stack’s part is a cheap, lazy look before an expensive ingest, since the source returns SourceRows without opening files and a coarse GeoSlice previews them. catalog_of stands for building an in-memory catalog from the rows (geocatalog.sources.source_row_to_gdf_row).

Case 18 · The reinforcement loop: verdicts become labels, labels retrain the model¶

TL;DR. Every analyst verdict and every dismissed alert becomes a training label, a frozen label set retrains the detector, and a promotion re-points the production tag.

Problem. The label sets $\mathcal{D}_1 \subset \mathcal{D}_2 \subset \ldots$ grow as verdicts arrive. We retrain on one frozen set (such as $\mathcal{D}_{12}$, the v12 of the code), compare the candidate with the production model on held-out scenes, and promote it only if it is better.

1
2
3
4
5
6
7
8
9
10
from pipekit_train import IterableDataset, TrainingLoop
from pipekit_experiment import LocalModelRegistry, TrainingArtifact

labels = frozen_label_set("v12")                     # 𝒟₁₂ from geoapps   user
ds = IterableDataset(source=lambda: chip_label_pairs(labels),
                     content_hash=labels.hash)       # today
candidate, _ = TrainingLoop(model_op=detector, dataset=ds).run()  # args sketched
ref = reg.store(candidate, weights=candidate.serialize_weights())   # today
if better(candidate, reg.load(reg.resolve_tag("plume-detector"))):  # user
    reg.tag(ref, "plume-detector", force=True)       # promotion             today

Notes. This is the loop that makes geoapps a model-reinforcement platform rather than a dashboard. The stack’s contribution is content-addressing at every step: the label set’s hash goes into the dataset, the dataset’s hash and the model’s hash go into a TrainingArtifact, and a tag names the production model. So a disputed detection can be traced to the exact labels its model was trained on. serialize_weights is pipekit_jax.JaxModelOp’s method, and other model wrappers bring their own.

Case 19 · Event kinds as pipelines¶

TL;DR. Support plumes, oil spills, floods and cyclones with one platform, where each kind is a pipeline and a few rules and never a branch in the core.

Problem. A kind is a contract made of 1) the operators that turn a scene into detections (retrieve → detect → quantify), 2) the marks a detection carries, 3) a link rule that chains detections into events, and 4) the entities it can be attributed to.

  EventKind  (one per phenomenon; a plugin)
    ├── operators: retrieve → detect → quantify   ← a pipekit Graph
    ├── marks:     what a detection measures         (flux ± σ, area, peak, …)
    ├── link rule: detections → events            (IoU carry-over, a tracker)
    └── entities:  what it can be attributed to
            │
  shared core  (kind-agnostic): scene · detection · event · coverage · label · run

Notes. Each kind’s operator chain (Table 4) is a pipekit Graph like the running example, so the platform stores, hashes, reruns and compares every kind the same way. For methane, the retrieve and quantify steps are the ones that the GeoModels boundaries move into radtrax and plumax. The link rule’s IoU carry-over matches the stack’s geocatalog.matchup.IouAtLeast, which pairs detections across overpasses. Rain and cyclone products are science-side xarray work rather than new app geometry, which is why xrtoolz carries them.

Case 20 · Reanalysis evaluation against in-situ (geocatalog → xrtoolz)¶

TL;DR. Score a gridded reanalysis against station observations, pairing each observation with the nearest grid time.

Problem. We have observations $y_j \in \mathbb{R}$ at stations $j = 1, \ldots, N_\text{obs}$, with locations $\vect{s}_j \in \mathbb{R}^2$ and times $t_j$, and a gridded product $x(\vect{s}, t) \in \mathbb{R}$. We pair each $y_j$ with $x(\vect{s}_j, t)$ at the nearest time within 3 hours and compute $\mathrm{RMSE} = \sqrt{\frac{1}{N_\text{obs}} \sum_j (x(\vect{s}_j, t_j) - y_j)^2}$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# Shapes: N_obs matched pairs; grid (time, lat, lon)
import pandas as pd
import pipekit_ops as pk
from geocatalog.matchup import Intersects, NearestInTime, matchup
from xrtoolz import geo, metrics

# pairs = {(j, t) : s_j ∈ cell(t), |t − t_j| ≤ 3 h}
pairs = list(matchup(stations, reanalysis, spatial=Intersects(),
                     temporal=NearestInTime(dt=pd.Timedelta("3h"))))   # N_obs rows  today
anom = pk.Sequential([
    geo.ValidateCoords(),                            # (time, lat, lon)     today
    # x̃ = x − x̄_doy
    geo.RemoveClimatology(clim),                     # (time, lat, lon)     today
])
# RMSE = √((1/N_obs) Σ_j (x(s_j, t_j) − y_j)²)
rmse = metrics.RMSE(variable="t2m", dims="obs")     # today
score = rmse(paired_dataset(pairs, anom))            # user: assemble the pairs

Notes. This is the xarray side that the original gallery lacked, and it uses the same geocatalog and pipekit as the methane cases, with xrtoolz in place of geotoolz. paired_dataset stands for assembling the matched pairs into one xr.Dataset with an obs dimension, which is user code today. The original’s hypothetical xrtoolz modules map onto real ones:

• detrend → xrtoolz.geo;

• calc → xrgrad;

• transforms → xrtoolz.transforms;

• ice → xrtoolz.ice, which is empty.

Which libraries each case needs¶

What the matrix (Table 5) says:

1. pipekit is in 19 of 20 cases. The exception is the explorer, which is pure catalog work, and it still belongs in a library of its own.

2. pipekit-experiment is the busiest sister package (10 cases), not pipekit-array (one optional use and one proposal), whereas the original gallery expected the reverse. Deployment needs hashing, YAML, registries and artifacts more than it needs array abstractions.

3. geocatalog is in 8 cases, geopatcher in 3. Both are used standalone, next to the domain libraries.

4. xrtoolz now has real cases (13, 19, 20), and they use the same framework and infrastructure as the methane cases. Swapping geotoolz for xrtoolz leaves the rest of each case intact, which is the strongest evidence that the framework layer is carrier-neutral.

5. The geoapps cases (14 to 19) lean on hashing and the registry. Versioned reruns, promotions and traceable labels are where a platform meets the stack.

Patterns the cases reinforce¶

pipekit owns composition; its sister packages own deployment infrastructure¶

Each concern still has one home (Table 6), so the original’s central claim holds, and the methane physics has its own owners in the GeoModels package boundaries. The difference is that pipekit core stayed small (about 3,200 lines), with the deployment infrastructure one package over. This answers the original’s second risk, which was that pipekit would accrete infrastructure code.

Three roles for georeader¶

Cases reach georeader in three ways (Table 7).

Parallelism is chosen per deployment¶

There is no automatic parallelism anywhere, by design, so each deployment picks its own primitive (Table 8).

YAML is the lingua franca, read-only for now¶

Cases 2, 4, 9, 12 and 13 load pipelines from YAML, and every library’s operators are reconstructible from _target_ paths because they all subclass pipekit_ops.Operator. So a YAML can name geotoolz.plume.IMEEstimate and xrtoolz.geo.RemoveClimatology side by side. Today the YAML is written by hand or generated with hydra-zen, and turning a built nested pipeline back into YAML is proposal A.

Where JAX helps¶

pipekit-jax has shipped, but not as the “JAX-traceable pipelines” library that the original gallery imagined (Report 3). It wraps an Equinox model as an operator with byte-identical weight round-trips, and it adds diffrax forward models for pipekit-cycle. For the cases here (Table 9):

The original gallery’s verdict stands: JAX changes how fast the inference-heavy cases run, but not how they are written. The pipeline itself is still eager Python.

Risks¶

• pipekit-array is behind. Only its first phase shipped. However, the cases need it in one place only (case 8, via proposal D), which is far less than planned. Fallback: pk.AssertCallable for the checks until phase C lands.

• Infrastructure code accretes. It did, but in pipekit-experiment rather than pipekit core, which kept core small. Fallback: keep each concern in its own module (adapters.hydra, registry, artifacts), as today.

• Glue reinvented in each consumer. CatalogPipeline-style glue (cases 2, 10, 12, 15) has no home, so every consumer writes its own loop. Recommended: proposal C, geocatalog.integrations.pipekit, mirroring geopatcher’s integration module. Alternative: a small pipekit-catalog package.

• YAML that only loads one way. A team that builds pipelines in Python cannot save them as YAML, which breaks cases 2, 5, 9 and 12 for nested pipelines. Fallback: keep the YAML as the source of truth until proposal A lands.

• Private API across packages. geotoolz and geopatcher’s pipekit bridge import pipekit_ops._base.operator helpers, and making them public is Report 16 §A.

Open questions¶

1. Does deployment infrastructure stay in pipekit-experiment? Hashing, YAML and artifacts are the busiest parts of the stack in these cases. Promoting a small pipeline_hash (proposal B) into pipekit core would make it reachable without the experiment package.

2. Is the geoapps job queue a stack concern? Case 15’s keyed, resumable steps look general (any batch job wants them). Today they are geoapps’ Postgres queue, so a small stack-level key-value “done” store might belong in pipekit-experiment.

3. Should event kinds be expressible as configs? If a kind’s operator chain is a YAML Graph, then the platform can add kinds without code. That depends on proposal A and on a YAML form of Graph, which does not exist.

4. How are human labels versioned? Case 18 assumes frozen label sets with a content hash, which pipekit-train’s datasets accept. So who owns the label store’s hashing: geoapps or pipekit-experiment?

Testing gates for the proposals¶

Each proposal merges only when its check passes (Table 10).

References¶

Sources¶

• The original use-case gallery, usecases.md.

• The geoapps design plan and repository (jejjohnson/geoapps), for cases 14 to 19.

• The pinned code in references/, and the runs described in each caption.