Use cases across the stack

The library shape, as it ships¶

Figure 1 is the stack the cases run on; its layering and ownership rules are set out in Report 16. One consequence 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 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 radiance $L \in \mathbb{R}^{B \times H \times W}$: $B$ spectral bands over an $H \times W$ pixel grid with a pixel area of $\Delta A$ (60 m pixels give $\Delta A = 3600 \mathrm{m}^2$). A methane plume absorbs in the shortwave infrared, so it shows up as a small, spectrally specific dip in $L$.

The enhancement. A matched filter against the background mean $\mu$, covariance $\Sigma$ and methane target signature $t$ estimates the column enhancement $\alpha(x)$ at each pixel $x$; its derivation and its limits are in radtrax · 1, stage 5.

The flux. Converting $\alpha$ to a mass column $\Omega(x)$ and masking the plume pixels $\mathcal{P}$ gives the integrated mass enhancement $\mathrm{IME} = \sum_{x \in \mathcal{P}} \Omega(x) \Delta A$, and the source rate $Q = U_{\mathrm{eff}} \cdot \mathrm{IME} / \ell$ for an effective wind speed $U_{\mathrm{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.

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 as pk
import geotoolz as gz

radiance = pk.Input("radiance")                          # L: (B, H, W)          today
# α(x) = (L(x) − μ)ᵀ Σ⁻¹ t / (tᵀ Σ⁻¹ t)
enh = gz.matched_filter.MatchedFilter(
    target=t, fit_on_call=True)(radiance)                # (B, H, W) → (H, W)   today
# Ω(x) = α(x) · ρ_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
# 𝒫 = {x : α(x) > τ_Otsu}, connected components ≥ 20 px
mask = gz.plume.PlumeMask(threshold="otsu", min_area_px=20)(enh)  # bool (H, W) today
# IME = Σ_{x∈𝒫} Ω(x) Δ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 most cases below pass around. Cases that need a single-input, single-output pipeline use a Sequential of the same operators.

What the cases need that does not ship yet¶

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

Exploration and batch processing¶

Case 1 · Notebook exploration (pipekit → geotoolz)¶

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

Problem. A Sentinel-2 scene $X \in \mathbb{R}^{B \times H \times W}$ with a QA60 band. Mask clouds, stretch each band, compute NDVI and see 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 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 + ε)
    gz.indices.NDVI(nir=7, red=3),                  # (B, H, W) → (1, H, W)  today
])
ndvi = viz(scene)                                   # GeoTensor          today

Library map. pipekit and geotoolz; the scene came from an earlier georeader call. Notes. 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, writing a COG per scene.

Problem. A catalog of $S$ scenes and a pipeline $f$ defined in YAML. Compute $f(X_s)$ for $s = 1 \dots S$ and write each result, tolerating individual failures.

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

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

# f(X_s) for s = 1 … S, 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), 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; pipekit.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 CI and production both fail loudly on bad data.

Problem. Before and after a correction step, 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 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)
    # f_valid = |{x : X(x) 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: they raise QCError, and Quarantine is the non-raising alternative. TOAToBOA from the original gallery does not exist; SimpleAtmosphericCorrection is the nearest.

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. Inputs $X_{\mathrm{opt}} \in \mathbb{R}^{B_o \times H \times W}$, $X_{\mathrm{sar}} \in \mathbb{R}^{B_s \times H' \times W'}$ and $Z \in \mathbb{R}^{H'' \times W''}$ on three grids. Align them to the optical grid, stack, 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 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; 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 Graph.summary() catches a missing input at construction. It does not exist: 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 $x_k \in \mathbb{R}^{B \times 256 \times 256}$, $k = 1 \dots K$. Inference splits a scene into $N_p$ patches at stride 224 (32 pixels of overlap), predicts $\hat y_p$ per patch, and stitches them 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_p inference patches
import pipekit 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, w) → (B, h, w)  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: y(x) = Σ_p w_p(x) ŷ_p(x) / Σ_p w_p(x)
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_p patches     today
    ApplyToChips(operator=preprocess | model),   # (B, 256²) → (1, 256²)  today
    MergePatches(aggregation=spatial.aggregation.OverlapAdd(),
                 domain=field.domain),               # N_p patches → (1, H, W)  today
])
prediction = infer(gp.RasterField(scene))           # today

Library map. geocatalog (chip sources), pipekit (composition), pipekit-train (the loop and dataset), geopatcher with its [pipekit] extra (patching), geotoolz (operators and ModelOp), georeader (the field). Notes. 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: 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 most uncertain hundred first.

Problem. For patches $p = 1 \dots N_p$, compute a per-patch score $s_p = \overline{\sigma_p}$, the mean predictive standard deviation, and select the $k = 100$ patches with the largest $s_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_p patches of 256 × 256; k = 100 selected
import heapq
import numpy as np
import pipekit 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
    # s_p = (1 / |Ω_p|) Σ_{x∈Ω_p} σ_p(x)
    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
)
# top-k of {s_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"). pipekit-array’s MeanScalar reduces over array axes (axis=); it has no field argument, so on a plain array pipekit_array.MeanScalar() works too. Scoring is the expensive part, so geopatcher.run.parallel_map is 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. 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, w)
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
        # h = 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, w)   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). 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 (checked). Proposal B would expose the hash without the throwaway registry. Heavy state lives in setup, never in predict. That is operational discipline, not 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 name → pipeline map. 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, w) → 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: point plume at a new hash, and every server that resolves the tag picks it up. One limit: 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, 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 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: loads_sandboxed refuses any class not registered. 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. 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 tile $(z, x, y)$, 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, w) 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; 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. 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, 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 records its class path.

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

TL;DR. Hand each scene to an orchestrator task that retries safely, writing 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 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 = H(pipeline) ‖ H(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; filter the catalog with query(catalog, slice_) for a date window first. The orchestrator owns retries and fan-out; 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). Report the config difference, and for a sample of $S$ scenes the per-scene $\mathrm{RMSE}_s = \lVert f(X_s) - f'(X_s) \rVert_2 / \sqrt{HW}$ and its distribution.

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: S 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_s = ‖f(X_s) − f′(X_s)‖₂ / √(HW),  s = 1 … S
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])            # (S,)                 user

Notes. A mean RMSE hides a bimodal regression, so look at the whole $\{\mathrm{RMSE}_s\}$ distribution, or better, the change in the quantity that matters ($\Delta Q_s$ 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; 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. Its design plan adds six shapes 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; 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 analysts validate the most urgent first; nothing unvalidated goes downstream.

Problem. Each predicted detection $d$ has a flux estimate $\hat Q_d$, an age $a_d$ and a flag $n_d \in \{0, 1\}$ (would validating it trigger a notification?). The priority combines them:

where $g$ increases with flux, $\tau$ is the notification window and $\lambda$ the decay after it. The code below computes (1) from the running example’s evidence.

1
2
3
4
5
6
# Shapes: D 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 + n_d) u(a_d)
priority = g(q_hat) * (1 + would_notify) * u(age)    # geoapps-core       user

Notes. The stack’s job ends at the evidence: an enhancement map, a mask, $\hat Q_d$ and $\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), and case 18 consumes it.

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 version tuple $v$ (pipeline hash, model hashes, input product versions). Each step’s output is keyed by $k = (\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: S 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
# k = (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; without them it is just a detection count. Reanalysis never overwrites: 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. Swap $U_{\mathrm{eff}}$ (ERA5) for $U'_{\mathrm{eff}}$ (GEOS-FP) with the mask held fixed. Since $Q = U_{\mathrm{eff}} \cdot \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: 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 window, list 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: 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. Label sets $\mathcal{L}_1 \subset \mathcal{L}_2 \subset \dots$ grow as verdicts arrive. Retrain on a frozen $\mathcal{L}_i$, 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. A disputed detection can then be traced to the exact labels its model was trained on. serialize_weights is pipekit_jax.JaxModelOp’s method; 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, never a branch in the core.

Problem. A kind is a contract: operators that turn a scene into detections (retrieve → detect → quantify), the marks a detection carries, a link rule that chains detections into events, and 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 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, not 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. Observations $y_j$ at stations $j = 1 \dots J$ and times $t_j$; a gridded product $x(s, t)$. Pair each $y_j$ with $x(s_j, t)$ at the nearest time within 3 hours, and compute $\mathrm{RMSE} = \sqrt{\frac{1}{J} \sum_j (x(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: J matched pairs; grid (time, lat, lon)
import pandas as pd
import pipekit 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"))))   # J rows  today
anom = pk.Sequential([
    geo.ValidateCoords(),                            # (time, lat, lon)     today
    # x′ = x − x̄_doy
    geo.RemoveClimatology(clim),                     # (time, lat, lon)     today
])
# RMSE = √((1/J) Σ_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 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. 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). 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; 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. That answers the original’s second risk, 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: each deployment picks its 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.Operator. 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; 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 the original gallery imagined (Report 3). It wraps an Equinox model as an operator with byte-identical weight round-trips, and 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, not how they are written. The pipeline itself is still eager Python.

Risks¶

• pipekit-array is behind. Only its first phase shipped. The cases need it in one place (case 8, via proposal D), 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._base.operator helpers; 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; 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, 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. 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.