Three kinds of function, three kinds of data

▦ Plain array

A block of numbers with a shape and nothing else:

• $N_\text{ax}$: the number of axes; $a = 1, \ldots, N_\text{ax}$ indexes them.
• $N_a$: the length of axis $a$.

Example The scene is $\vect{x} \in \mathbb{R}^{N_\text{band} \times N_\text{row} \times N_\text{col}}$, so $N_\text{ax} = 3$. Which axis is the band is a convention the code has to remember. EMIT's radiance file stores (y, x, band) and geoarray (band, y, x); a mean over axis 0 meant as the mean spectrum averages down the rows instead, and runs without an error.

In code a NumPy or JAX array. Its axis order lives in a jaxtyping annotation and in each einx pattern, so the convention is written down at every call, but nothing checks it against the data: a scene read as (y, x, band) passes the same call and averages over the wrong axes.

Reading the sketches. Each block opens with a # Shapes: line naming its dimensions; the column on the right tracks shapes, (B, R, C) → (B,), in those letters. A call into the stack has a marker above it: # today (pkg) ships now, # proposed (pkg) does not yet, and # user is a function or class this page defines. The equation it computes sits above each line.

import einx
import jax.numpy as jnp
import numpy as np
from jaxtyping import Array, Float

# Shapes: B = N_band = 285 bands; R × C = N_row × N_col = 256 × 256 pixels


def mean_spectrum(radiance: Float[Array, "band y x"]) -> Float[Array, " band"]:
    """Mean radiance spectrum of a scene, averaged over its pixels.

        x̄_b = (1 / R C) ∑_ij x_bij

    Algorithm:
        1. Average over the axes the pattern names y and x.

    Args:
        radiance: x, the scene, with the band axis first.

    Returns:
        x̄, one value per band. Nothing checks that the band axis really is
        first: a (y, x, band) array gives a wrong answer and no error.
    """
    # x̄_b = (1 / R C) ∑_ij x_bij; the axis names live in this pattern only
    spectrum = einx.mean("band [y x]", radiance)        # (B, R, C) → (B,)
    return spectrum


# x ∈ ℝ^(B × R × C): a crop of an EMIT scene, on the host
radiance = np.zeros((285, 256, 256), dtype=np.float32)  # → (B, R, C)
# the same numbers on the device
scene = jnp.asarray(radiance)                           # (B, R, C) → (B, R, C)
# user
# x̄
spectrum = mean_spectrum(scene)                         # (B, R, C) → (B,)

▤ Named array

The same numbers with a name $d_a$ and a coordinate vector $\vect{c}_a$ for each axis:

• $\vect{d}$: one name per axis, all distinct.
• $\vect{c}_a \in \mathbb{R}^{N_a}$: the coordinates along axis $a$ (wavelengths, positions, dates).

An operation addresses an axis by name, and the array finds its position. Take a reduction $r_a$ (a mean, a sum, a maximum) over axis $a$ of a plain array: it maps $\mathbb{R}^{N_1 \times \ldots \times N_{N_\text{ax}}}$ to the same space with the factor $N_a$ left out. On a named array the same reduction is asked for by a name $d$. The array looks up its position, $a = \op{pos}_{\vect{d}}(d)$, reduces there, and removes that axis's name and coordinates, since the output no longer has that axis:

• $\op{pos}_{\vect{d}}(d)$: the position of the name $d$ in $\vect{d}$; an error if $d$ is not there.
• $\vect{d}_{-a}$, $\mathcal{C}_{-a}$: the names and coordinates with entry $a$ taken out.

Example $\vect{d} = (\texttt{band}, \texttt{y}, \texttt{x})$, with $\vect{c}_1$ the band wavelengths. The mean spectrum reduces by name twice, over $\texttt{y}$ at $\op{pos}_{\vect{d}}(\texttt{y}) = 2$ and then over $\texttt{x}$, and returns names $(\texttt{band})$ with the wavelengths still attached, whatever order the reader used.

In code pipekit-array's DimArray: the data, one name per axis, and static labels on the dims that have them (here the wavelengths), with no other coordinates. It is a JAX pytree whose only leaf is the data, so new values never recompile. An xarray DataArray is the same object with full coordinates, and pipekit's dim functions (reduce, sel, isel, transpose) run on both through their adapters.

import einx
import numpy as np
import pipekit_array as pa

# Shapes: B = N_band = 285 bands; R × C = N_row × N_col = 256 × 256 pixels

radiance = np.zeros((285, 256, 256), dtype=np.float32)  # → (B, R, C)
# c_1: the centre wavelength of each band, nm
wavelengths = np.linspace(381.0, 2493.0, 285)           # → (B,)
rounded = np.round(wavelengths, 1)                      # (B,) → (B,)
WAVELENGTHS_NM = tuple(rounded.tolist())                # (B,) → B labels

# today (pipekit-array)
# x_named = (x, d, C): one name per axis, the band dim labelled by wavelength
scene = pa.DimArray(radiance, ("band", "y", "x"), {"band": WAVELENGTHS_NM})
# today (pipekit-array)
# r_y, then r_x: the mean spectrum, found by name
spectrum = pa.reduce(scene, "mean", ("y", "x"))         # (B, R, C) → (B,)
# the file's own order, (y, x, band)
as_stored = einx.id("band y x -> y x band", radiance)   # (B, R, C) → (R, C, B)
# today (pipekit-array)
# the same numbers, named in that order
stored = pa.DimArray(as_stored, ("y", "x", "band"), {"band": WAVELENGTHS_NM})
# today (pipekit-array)
# r_y, then r_x: the same call gives the same answer
same_spectrum = pa.reduce(stored, "mean", ("y", "x"))   # (R, C, B) → (B,)

◆ Carrier

A named array plus its georeference $\Gamma$:

• $\op{crs}$: the coordinate reference system (an EPSG code or WKT).
• $\mat{A} \in \mathbb{R}^{2 \times 3}$: the affine transform from pixel to map coordinates.
• $x_\text{fill} \in \mathbb{R}$: the fill value, NaN allowed.
• $\mathcal{M}$: the attributes (units, sensor, the sun and view zenith angles, processing level). Band names and wavelengths are coordinates of the band axis, in $\mathcal{C}$.

The transform says where each pixel is. The top-left corner of the pixel in row $i$ and column $j$ sits at map coordinates $\vect{s}_{ij} \in \mathbb{R}^2$, and its centre at $(j + 1/2, i + 1/2)$ in pixel units:

The fill value says which pixels are missing. A pixel is valid when every band there is finite and differs from the fill, which gives the valid mask $\vect{m} \in \{0, 1\}^{N_\text{row} \times N_\text{col}}$, with $b$ indexing bands:

For a swath or a curvilinear grid, as EMIT's radiance is before orthorectification, $\mat{A}$ is replaced by per-pixel coordinates $\mat{S}_\text{lon}, \mat{S}_\text{lat} \in \mathbb{R}^{N_\text{row} \times N_\text{col}}$; the rest is the same.

Example The scene is a carrier on a 60 m grid in UTM zone 13N with $x_\text{fill} = -9999$, the band names and centre wavelengths in $\mathcal{C}$, and the sun and view zenith angles in $\mathcal{M}$. The well pad at map coordinates $\vect{s}$ is found at pixel $(i, j)$ by solving $\vect{s} = \mat{A} (j, i, 1)^\top$ for $(i, j)$, which is to_pixel. Slice the bare array instead, np.asarray(scene)[..., 100:164, 100:164], and $\mat{A}$ still describes the full scene: the window claims to start at the scene's corner, and the plume is attributed to a facility several kilometres away.

In code geotoolz-array's GeoArray: the data (NumPy or JAX), its dims ending in ("y", "x"), coordinates for the other dims (SpectralAxis, TimeAxis), the transform, CRS, fill value and frozen attrs. Slicing, crop and read_window move the transform with the pixels; a reduction over y and x leaves no grid, so it returns a named array. GeoSwath carries per-pixel latitude and longitude instead of a transform, and GeoPoints a set of points; xrops plays the same role for a DataArray whose coordinates and CRS mean geography.

import numpy as np
import pipekit_array as pa
from geoarray import GeoArray, crop, to_pixel
from geoarray.coords import SpectralAxis

# Shapes: B = N_band = 285 bands; R × C = N_row × N_col = 256 × 256 pixels

radiance = np.zeros((285, 256, 256), dtype=np.float32)  # → (B, R, C)
# the centre wavelength of each band, nm, and a name for each
wavelengths = np.linspace(381.0, 2493.0, 285)           # → (B,)
rounded = np.round(wavelengths, 1)                      # (B,) → (B,)
WAVELENGTHS_NM = tuple(rounded.tolist())                # (B,) → B labels
BAND_NAMES = tuple(f"{wavelength:.1f}nm" for wavelength in WAVELENGTHS_NM)

# today (geoarray)
# C: the band names with their centre wavelengths
bands = SpectralAxis(names=BAND_NAMES, wavelength_nm=WAVELENGTHS_NM)
# today (geoarray)
# x_geo = (x_named, Γ), Γ = (crs, A, x_fill, M): 60 m in UTM zone 13N
scene: GeoArray["band y x"] = GeoArray(                 # (B, R, C) → (B, R, C)
    radiance,
    (60.0, 0.0, 630000.0, 0.0, -60.0, 3550000.0),
    "EPSG:32613",
    fill_value=-9999.0,
    dims="band y x",
    coords=(bands,),
    attrs={"solar_zenith_deg": 28.4, "view_zenith_deg": 6.1},
)
# today (geoarray)
# m_ij = ∏_b 1[x_bij ≠ x_fill] 1[x_bij is finite]
valid = scene.valid_mask()                              # (B, R, C) → (R, C)
# s: the well pad in map coordinates, m
well_pad_xy = np.array([[637530.0, 3542210.0]])         # → (1, 2)
# today (geoarray)
# s = A (j + ½, i + ½, 1)ᵀ, solved for the fractional (i, j)
well_pad_pixel = to_pixel(scene.transform, well_pad_xy)  # (1, 2) → (1, 2)
# today (geoarray)
# a 64 × 64 window from row 100, column 100; A moves with the pixels
window = crop(scene, (100, 100), (64, 64))              # (B, R, C) → (B, 64, 64)
# today (pipekit-array)
# r_y, then r_x: no grid is left, so a named array comes back
spectrum = pa.reduce(scene, "mean", ("y", "x"))         # (B, R, C) → (B,)

ƒ Function

A pure function of arrays, with its parameters passed in on every call:

• $\vect{x}$: the input array; $\vect{y}$: the output array.
• $\vect{u}$: what $f$ needs to know about where the data sits, as plain arrays or numbers (pixel areas, a transform, an air-mass factor, absorption at the band wavelengths).
• $\gvect{\theta} \in \mathbb{R}^{D_\theta}$: parameters, constant over the call.

$f$ never sees a name, a CRS or a fill value. That is what lets it run under jit, on any carrier, in any library below GeoStack.

Example The matched filter, the first step of the running example. With $n = 1, \ldots, N_\text{row} N_\text{col}$ indexing pixels and $\vect{x}_n \in \mathbb{R}^{N_\text{band}}$ the spectrum at pixel $n$, it estimates the methane enhancement $y_n$, in ppm m, that best explains how $\vect{x}_n$ departs from the background. Here the background is given, as the parameters $\gvect{\theta} = (\gvect{\mu}, \gmat{\Sigma})$:

The target signature $\vect{t}$ is the radiance change for one ppm m of methane. Linearising the Beer–Lambert law, $x_{n,b} \approx \mu_b (1 - M k_b y_n)$ for band $b$, gives

• $\gvect{\mu} \in \mathbb{R}^{N_\text{band}}$, $\gmat{\Sigma} \in \mathbb{R}^{N_\text{band} \times N_\text{band}}$: the background mean and covariance.
• $\vect{k} \in \mathbb{R}^{N_\text{band}}$: the methane absorption at each band's wavelength, per ppm m along a vertical path.
• $M \in \mathbb{R}$: the geometric air-mass factor, one over the cosine of the solar zenith angle plus one over the cosine of the view zenith angle.
• $\vect{u} = (M, \vect{k})$: both read off the data, the angles from the attributes and the absorption at the band coordinate's wavelengths.
• $\vect{t} \in \mathbb{R}^{N_\text{band}}$: the target signature; $\odot$ multiplies element by element.

With the background given, the estimate is the generalised least-squares fit of a simple model, a pixel's spectrum as the background plus $y_n$ ppm m of methane along $\vect{t}$, with Gaussian clutter $\gvect{\varepsilon}_n$; it is also the maximum-likelihood estimate:

Setting the derivative in $y$ to zero gives the closed form above; its standard error is $\sigma_y = (\vect{t}^\top \gmat{\Sigma}^{-1} \vect{t})^{-1/2}$, the same at every pixel.

The output keeps a singleton band axis, $\vect{y} \in \mathbb{R}^{1 \times N_\text{row} \times N_\text{col}}$.

import einx
import lineax as lx
from jaxtyping import Array, Float

# Shapes: B = N_band bands; R × C = N_row × N_col pixels


def matched_filter(
    radiance: Float[Array, "band y x"],
    unit_absorption: Float[Array, " band"],
    air_mass: float,
    *,
    mean: Float[Array, " band"],
    covariance: lx.AbstractLinearOperator,
) -> Float[Array, "1 y x"]:
    """Methane enhancement of every pixel by a matched filter, background given.

    Model: x_n = μ + y_n t + ε_n, ε_n ~ 𝒩(0, Σ), with θ = (μ, Σ) known. The
    estimate is the generalised least-squares fit of y_n, which is also its
    maximum-likelihood estimate:

        y_n = argmin_y (x_n − μ − y t)ᵀ Σ⁻¹ (x_n − μ − y t)
            = tᵀ Σ⁻¹ (x_n − μ) / (tᵀ Σ⁻¹ t),  t = −M (k ⊙ μ)

    Algorithm:
        1. t = −M (k ⊙ μ): the linearised Beer–Lambert response to 1 ppm m.
        2. s = Σ⁻¹ t, by a solve; a function repeats it on every call.
        3. y_n = sᵀ (x_n − μ) / (tᵀ s) at every pixel.

    Args:
        radiance: x, the scene's spectra.
        unit_absorption: k, methane absorption per ppm m at each band; part of u.
        air_mass: M, the geometric air-mass factor; part of u.
        mean: μ, the background mean; part of θ.
        covariance: Σ, the background covariance as a PSD operator; part of θ.

    Returns:
        y in ppm m, with a singleton band axis; its standard error is
        (tᵀ Σ⁻¹ t)^(−1/2) at every pixel.
    """
    # k ⊙ μ
    absorbed = unit_absorption * mean                   # (B,), (B,) → (B,)
    # t = −M (k ⊙ μ): the radiance change for 1 ppm m of methane
    target = -air_mass * absorbed                       # (B,) → (B,)
    # s = Σ⁻¹ t
    solution = lx.linear_solve(covariance, target)      # (B, B) PSD op, (B,)
    filter_direction = solution.value                   # → (B,)
    # tᵀ s
    normalisation = einx.dot("band, band ->", target, filter_direction)  # → ()
    # r_n = x_n − μ
    residual = einx.subtract(                           # (B, R, C), (B,) → (B, R, C)
        "band y x, band -> band y x", radiance, mean
    )
    # sᵀ r_n
    projection = einx.dot(                              # (B, R, C), (B,) → (R, C)
        "band y x, band -> y x", residual, filter_direction
    )
    # y_n = sᵀ r_n / (tᵀ s)
    enhancement = projection / normalisation            # (R, C) → (R, C)
    # a singleton band axis
    return einx.id("y x -> 1 y x", enhancement)         # (R, C) → (1, R, C)

▣ Module

A function with state. A setup step builds the state $\gvect{\psi}$ once, from data $\mat{X}$ and what it needs to know about them, $\vect{u}$, and may do work that cannot be traced (a scikit-learn estimator, a k-NN graph, a factorisation). The call is then pure and jittable:

• $\mat{X}$: the data setup sees, one row per sample; for the matched filter, the spectra of the valid pixels, $\mat{X} \in \mathbb{R}^{N_\text{valid} \times N_\text{band}}$.
• $\gvect{\psi}$: the state, a pytree of arrays held as the Module's fields next to $\gvect{\theta}$.

Example The matched filter as it is used: the plume covers a small fraction of the pixels, so the scene is its own background. Setup estimates $\gvect{\mu}$ and a shrunk $\gmat{\Sigma}$ from $\mat{X}$ (Ledoit–Wolf, in scikit-learn), forms $\vect{t}$, and solves once for the filter direction $\vect{s} = \gmat{\Sigma}^{-1} \vect{t}$; the state is $\gvect{\psi} = (\gvect{\mu}, \vect{s}, \vect{t}^\top \vect{s})$. Each call is then two contractions:

What each method solves. Both rest on one model of a pixel's spectrum, the background plus $y_n$ ppm m of methane along $\vect{t}$:

fit solves for the background and the filter, once per scene and off the trace: a covariance shrunk towards a scaled identity by the amount that minimises its expected error, and the direction that maximises the signal-to-clutter ratio,

The call solves for each pixel's enhancement with that state fixed. It is the generalised least-squares fit of the model, which is also the maximum-likelihood estimate under the Gaussian background:

• $\gvect{\varepsilon}_n \in \mathbb{R}^{N_\text{band}}$: the background's departure at pixel $n$, the clutter the filter suppresses.
• $\mat{S} \in \mathbb{R}^{N_\text{band} \times N_\text{band}}$: the sample covariance of $\mat{X}$; $\nu = \op{tr}(\mat{S}) / N_\text{band}$; $\rho \in [0, 1]$: the Ledoit–Wolf shrinkage, the value that minimises $\mathbb{E} \Vert \hat{\gmat{\Sigma}} - \gmat{\Sigma} \Vert_F^2$.
• The argmin has the closed form above, $\vect{s}^\top (\vect{x}_n - \gvect{\mu}) / (\vect{t}^\top \vect{s})$, with standard error $\sigma_y = (\vect{t}^\top \vect{s})^{-1/2}$ at every pixel.

The docstrings below carry the same split: each method states its objective, its steps and its arguments.

import einx
import equinox as eqx
import jax.numpy as jnp
import numpy as np
import scipy.linalg
import sklearn.covariance
from jaxtyping import Array, Float

# Shapes: B = N_band bands; R × C = N_row × N_col pixels; P = N_valid valid pixels


class MatchedFilter(eqx.Module):
    """Matched filter whose background is fitted once, from the scene's own pixels.

    Model: a pixel's spectrum is the background plus y ppm m of methane along t,

        x_n = μ + y_n t + ε_n,  ε_n ~ 𝒩(0, Σ)

    `fit` estimates the background and the filter that model implies; the call
    estimates y_n at every pixel with them. The state is ψ = (μ, s, tᵀ s).
    """

    mean: Float[Array, " band"]
    filter_direction: Float[Array, " band"]
    normalisation: Float[Array, ""]

    @classmethod
    def fit(
        cls,
        spectra: Float[np.ndarray, "pixel band"],
        unit_absorption: Float[np.ndarray, " band"],
        air_mass: float,
    ) -> "MatchedFilter":
        """Fit the background and the filter direction, once per scene.

        Two objectives, solved in turn. The background covariance is shrunk
        towards a scaled identity by the amount that minimises its expected
        squared error,

            Σ̂ = (1 − ρ) S + ρ ν I,  ρ = argmin_ρ 𝔼 ‖Σ̂ − Σ‖²_F

        and the filter direction maximises the signal-to-clutter ratio,

            s = argmax_s (sᵀ t)² / (sᵀ Σ̂ s) ∝ Σ̂⁻¹ t,  t = −M (k ⊙ μ̂)

        Algorithm:
            1. μ̂, Σ̂: the sample mean and the Ledoit–Wolf shrunk covariance of X.
            2. t = −M (k ⊙ μ̂): the radiance change for 1 ppm m of methane.
            3. s = Σ̂⁻¹ t, one solve, and the normalisation tᵀ s.

        Args:
            spectra: X, the spectra of the valid pixels, one row each; the plume
                is a small fraction of them, so the scene is its own background.
            unit_absorption: k, methane absorption per ppm m at each band.
            air_mass: M, the geometric air-mass factor.

        Returns:
            The filter with its state ψ = (μ̂, s, tᵀ s); S is the sample
            covariance, ν = tr(S) / B and ρ ∈ [0, 1] the shrinkage.
        """
        # --- setup: host-side NumPy, SciPy and scikit-learn, runs once ---
        # μ̂, Σ̂ = (1 − ρ) S + ρ ν I
        background = sklearn.covariance.LedoitWolf().fit(spectra)  # (P, B) → fitted
        mean = background.location_                     # → (B,)
        covariance = background.covariance_             # → (B, B)
        # k ⊙ μ̂
        absorbed = unit_absorption * mean               # (B,), (B,) → (B,)
        # t = −M (k ⊙ μ̂)
        target = -air_mass * absorbed                   # (B,) → (B,)
        # s = Σ̂⁻¹ t, solved once
        filter_direction = scipy.linalg.solve(          # (B, B), (B,) → (B,)
            covariance, target, assume_a="pos"
        )
        # tᵀ s
        normalisation = einx.dot("band, band ->", target, filter_direction)  # → ()
        # ψ = (μ̂, s, tᵀ s), as device arrays
        return cls(
            mean=jnp.asarray(mean),
            filter_direction=jnp.asarray(filter_direction),
            normalisation=jnp.asarray(normalisation),
        )

    def __call__(self, radiance: Float[Array, "band y x"]) -> Float[Array, "1 y x"]:
        """Estimate the methane enhancement at every pixel, with ψ fixed.

        The generalised least-squares fit of the model to each pixel, which is
        the maximum-likelihood estimate under the Gaussian background:

            y_n = argmin_y (x_n − μ − y t)ᵀ Σ⁻¹ (x_n − μ − y t)
                = sᵀ (x_n − μ) / (tᵀ s)

        with standard error σ_y = (tᵀ s)^(−1/2), the same at every pixel.

        Algorithm:
            1. r_n = x_n − μ: the departure from the background.
            2. sᵀ r_n: its projection on the filter direction.
            3. y_n = sᵀ r_n / (tᵀ s): scaled to ppm m.

        Args:
            radiance: x, the scene's spectra, any number of pixels.

        Returns:
            y in ppm m, with a singleton band axis. Pure and jittable: no fit.
        """
        # --- every call: two contractions, pure and jittable ---
        # r_n = x_n − μ
        residual = einx.subtract(                       # (B, R, C), (B,) → (B, R, C)
            "band y x, band -> band y x", radiance, self.mean
        )
        # sᵀ r_n
        projection = einx.dot(                          # (B, R, C), (B,) → (R, C)
            "band y x, band -> y x", residual, self.filter_direction
        )
        # y_n = sᵀ r_n / (tᵀ s)
        enhancement = projection / self.normalisation   # (R, C) → (R, C)
        # a singleton band axis
        return einx.id("y x -> 1 y x", enhancement)     # (R, C) → (1, R, C)

→ Operator

A pipeline step that lifts $f$ onto carriers. The Operator $h[f, \gvect{\theta}]$ is built from the function and its parameters; it is called on a carrier with its config $\alpha$, a JSON record (which bands, which dimension, the fill policy). It unwraps the carrier with $\pi_\alpha$, calls $f$ on plain arrays, and rewraps with $\iota_\alpha$:

Read the first line from the inside out: $\pi_\alpha$ unwraps, $f$ computes, $\iota_\alpha$ rewraps. The unwrap hands down the array and what $f$ needs to know about it, $\vect{u}$, with every name in $\alpha$ already resolved against the carrier's names and coordinates. $\gvect{\theta}$ passes through untouched; $\alpha$ never reaches $f$. The rewrap builds the output carrier:

• $\alpha$: the config, names and settings only, so it round-trips through JSON.
• $\pi_\alpha$, $\iota_\alpha$: the unwrap and the rewrap, both read $\alpha$; $\iota_\alpha$ also takes the input carrier, whose georeference the output inherits or regrids.
• $\vect{d}_y$, $\mathcal{C}_y$, $\Gamma_y$: the output's names, coordinates and georeference.
• $\op{regrid}_\alpha$: the input grid to the output grid, for a reprojection or a resampling.

Operators compose: $h[f_2, \gvect{\theta}_2](h[f_1, \gvect{\theta}_1](\vect{x}_\text{geo}; \alpha_1); \alpha_2)$ is what pipekit writes step1 | step2, each step holding its own $\gvect{\theta}$ and $\alpha$.

Example The settings are the band dim and a spectral window, $\alpha = \{\text{dim}: \texttt{band}, \text{window}: [2122, 2488] \text{ nm}\}$. $\pi_\alpha$ resolves them against the carrier: it keeps the bands whose wavelengths in $\mathcal{C}$ fall in the window, computes $\vect{u} = (M, \vect{k})$ from the angles in $\mathcal{M}$ and those wavelengths, and stacks $\mat{X}$ from the pixels where $[\vect{m}]_{ij} = 1$. The Module's setup runs on $\mat{X}$, its call on the window. $\iota_\alpha$ returns a carrier on the same grid, $\Gamma_y = \Gamma$, with one band named ch4 and NaN where $[\vect{m}]_{ij} = 0$. geoops ships this step as geoops.matched_filter.MatchedFilter, which fits its background on the first call.

import einx
import numpy as np
import pipekit
from geoarray import GeoArray
from geoops.carrier import wrap_like
from radtrax.atmosphere import air_mass_factor

# Shapes: B = N_band bands; W window bands; R × C = N_row × N_col pixels; P valid

# MatchedFilter is the Module above; ch4_unit_absorption is proposed for radtrax


class MethaneEnhancement(pipekit_ops.Operator):
    """CH₄ enhancement on a carrier: the MatchedFilter Module lifted onto GeoArrays.

        h[f, ψ](x_geo; α) = ι_α(f(π_α(x_geo); ψ), x_geo)

    α = (dim, window_nm) says what to read; π_α resolves it against the
    carrier's coordinates C and attrs 𝓜, f is the Module, ι_α rewraps.
    """

    def __init__(
        self, *, dim: str = "band", window_nm: tuple[float, float] = (2122.0, 2488.0)
    ) -> None:
        """Hold α: names and settings only, so the config round-trips through JSON.

        Args:
            dim: the name of the band dimension.
            window_nm: the spectral window, nm, where methane absorbs strongly.
        """
        self.dim = dim
        self.window_nm = window_nm

    def _apply(self, scene: GeoArray["band y x"]) -> GeoArray["1 y x"]:
        """Unwrap, fit the filter, filter, rewrap.

            π_α(x_geo) = (x_W, u, X),  u = (M, k)
            ψ = setup(X, u),  y = f(x_W; ψ)
            y_geo = ι_α(y, x_geo),  Γ_y = Γ

        Algorithm:
            1. x_W: the bands whose centre wavelengths fall inside the window.
            2. u = (M, k): M from the sun and view angles in 𝓜, k at the
               window's wavelengths.
            3. X: the spectra of the pixels where the valid mask m is 1.
            4. Fit the Module on X and u, then call it on x_W.
            5. Rewrap on the same grid: one band named ch4, NaN where m = 0.

        Args:
            scene: x_geo, the radiance, with band wavelengths in C and
                solar_zenith_deg, view_zenith_deg in 𝓜.

        Returns:
            y_geo, the enhancement in ppm m on the scene's grid.
        """
        # --- π_α: the names in α, resolved against C and M ---
        # today (geoarray)
        # the band coordinate: names and centre wavelengths
        bands = scene.coord(self.dim)
        lower, upper = self.window_nm
        # the bands whose wavelengths fall inside the window
        window = [
            name
            for name, wavelength in zip(bands.names, bands.wavelength_nm)
            if lower <= wavelength <= upper
        ]
        # today (geoarray)
        swir = scene.sel({self.dim: window})            # (B, R, C) → (W, R, C)
        # today (radtrax)
        # M = 1 / cos ϑ_sun + 1 / cos ϑ_view, from the angles in the attrs 𝓜
        air_mass = air_mass_factor(                     # → ()
            swir.attrs["solar_zenith_deg"], swir.attrs["view_zenith_deg"]
        )
        # today (geoarray)
        wavelengths = swir.coord(self.dim).wavelength_nm  # → W labels
        # proposed (radtrax)
        # k: methane absorption per ppm m at each band's wavelength
        unit_absorption = ch4_unit_absorption(wavelengths)  # → (W,)
        # today (geoarray)
        # m: valid where every band is finite and differs from the fill
        valid = swir.valid_mask()                       # (W, R, C) → (R, C)
        # X: the spectra of the valid pixels, one row each
        swir_values = np.asarray(swir)                  # (W, R, C) → (W, R, C)
        pixels = einx.id("band y x -> (y x) band", swir_values)  # → (R C, W)
        is_valid = einx.id("y x -> (y x)", valid)       # (R, C) → (R C,)
        spectra = pixels[is_valid]                      # (R C, W) → (P, W)
        # --- the Module: setup once, then a pure call on plain arrays ---
        # user
        # ψ = setup(X, u)
        matched_filter = MatchedFilter.fit(spectra, unit_absorption, air_mass)
        # user
        # y = f(x; ψ)
        enhancement = matched_filter(swir.data)         # (W, R, C) → (1, R, C)
        # --- ι_α: same grid, one band in ppm m, NaN where m = 0 ---
        # y_n = NaN where m_n = 0
        masked = np.where(valid, enhancement, np.nan)   # (1, R, C) → (1, R, C)
        # today (geoops)
        return wrap_like(swir, masked, band_names=("ch4",), fill_value=np.nan)

The sketches follow the house style: shapes as jaxtyping types, the equation above each line, einx for reshapes and contractions. All six run as written (checked on 6 Nov 2026): the data-axis examples against geotoolz-array and pipekit-array, and the function, the Module and the Operator on a synthetic scene with lineax and scikit-learn. In the Operator, air_mass_factor is radtrax's (radtrax.atmosphere) and ch4_unit_absorption stands for a lookup table of methane absorption (radtrax, or plumax's hapi_lut); geoops ships the matching Operators, matched_filter.MatchedFilter, plume.PlumeMask, plume.ColumnToMass and plume.IMEEstimate.

⓪ The data

The scene and the wind are carriers, each on its own grid: the radiance on 60 m in UTM zone 13N with the sun and view angles in its attributes $\mathcal{M}$, the wind on ERA5's 0.25° grid in longitude and latitude. Nothing in the rungs above builds a grid or reads a file.

import numpy as np
from geoarray import GeoArray
from geoarray.coords import SpectralAxis

# Shapes: B = 285 bands; R × C = 256 × 256 pixels; G × H = 8 × 8 wind cells

# x: radiance from the orthorectified EMIT file
radiance = np.zeros((285, 256, 256), dtype=np.float32)  # → (B, R, C)
# v: ERA5 10 m wind speed around the field, one band on a 0.25° grid
wind_speed = np.full((1, 8, 8), 3.5)                    # → (1, G, H)
# the centre wavelength of each band, nm, and a name for each
wavelengths = np.linspace(381.0, 2493.0, 285)           # → (B,)
rounded = np.round(wavelengths, 1)                      # (B,) → (B,)
WAVELENGTHS_NM = tuple(rounded.tolist())                # (B,) → B labels
BAND_NAMES = tuple(f"{wavelength:.1f}nm" for wavelength in WAVELENGTHS_NM)

# today (geoarray)
# C: the band names with their centre wavelengths
bands = SpectralAxis(names=BAND_NAMES, wavelength_nm=WAVELENGTHS_NM)
# today (geoarray)
# x_geo: a 60 m grid in UTM zone 13N, with the sun and view angles in 𝓜
scene: GeoArray["band y x"] = GeoArray(                 # (B, R, C) → (B, R, C)
    radiance,
    (60.0, 0.0, 630000.0, 0.0, -60.0, 3550000.0),
    "EPSG:32613",
    fill_value=-9999.0,
    dims="band y x",
    coords=(bands,),
    attrs={"solar_zenith_deg": 28.4, "view_zenith_deg": 6.1},
)
# today (geoarray)
# v_geo: the wind on its own grid, in longitude and latitude
wind_field: GeoArray["band y x"] = GeoArray(            # (1, G, H) → (1, G, H)
    wind_speed, (0.25, 0.0, -105.0, 0.0, -0.25, 33.0), "EPSG:4326", dims="band y x"
)

① Primitives

Functions on arrays that know nothing of methane plumes: the target signature, the area of a pixel of an affine grid, and the large connected components of a mask.

• $a \in \mathbb{R}$: the pixel area in m², one value on an affine grid; $\mat{A}_{2 \times 2}$ is the linear part of $\mat{A}$.

The third primitive keeps only the large connected pieces of a candidate mask $\vect{c}$, with $K(n)$ the 4-connected component of $\vect{c}$ that holds pixel $n$:

import einx
import numpy as np
import scipy.ndimage
from jaxtyping import Array, Bool, Float

# Shapes: W = window bands; R × C = N_row × N_col pixels; K labelled components


def target_signature(
    mean: Float[Array, " band"], unit_absorption: Float[Array, " band"], air_mass: float
) -> Float[Array, " band"]:
    """The radiance change for 1 ppm m of methane: the matched filter's target.

    Linearising the Beer–Lambert law around the background,
    x_b ≈ μ_b exp(−M k_b y) ≈ μ_b (1 − M k_b y), so ∂x/∂y at y = 0 is

        t = −M (k ⊙ μ)

    Args:
        mean: μ, the background mean spectrum.
        unit_absorption: k, absorption per ppm m along a vertical path.
        air_mass: M, the geometric air-mass factor.

    Returns:
        t, in radiance units per ppm m; negative where methane absorbs.
    """
    # k ⊙ μ
    absorbed = unit_absorption * mean                   # (W,), (W,) → (W,)
    # t = −M (k ⊙ μ)
    target = -air_mass * absorbed                       # (W,) → (W,)
    return target


def pixel_area(transform: Float[np.ndarray, " 6"]) -> float:
    """The area of every pixel of an affine grid.

        a = |det A₂ₓ₂| = |a₀ a₄ − a₁ a₃|

    for A = (a₀, a₁, a₂, a₃, a₄, a₅), the transform in (a, b, c, d, e, f) order.

    Args:
        transform: A, six numbers, in map units per pixel.

    Returns:
        a, in squared map units (m² on a projected grid), the same at every pixel.
    """
    # det A₂ₓ₂ = a₀ a₄ − a₁ a₃
    determinant = transform[0] * transform[4] - transform[1] * transform[3]  # → ()
    # a = |det A₂ₓ₂|
    return abs(determinant)


def large_components(
    mask: Bool[np.ndarray, "y x"], min_pixels: int
) -> Bool[np.ndarray, "y x"]:
    """The connected components of a mask that have at least N_min pixels.

        b_n = c_n 1[|K(n)| ≥ N_min]

    with K(n) the 4-connected component of c that holds pixel n.

    Algorithm:
        1. Label the components of c 1, …, K, and 0 off the mask.
        2. Count the pixels of each label.
        3. Keep the labels with at least N_min pixels; never label 0.

    Args:
        mask: c, the candidate pixels, above a threshold.
        min_pixels: N_min.

    Returns:
        b, the mask with its small components removed.
    """
    # the components of c, labelled 1, …, K; 0 off the mask
    labels, count = scipy.ndimage.label(mask)           # (R, C) → (R, C), K
    # |component k|, in pixels
    flat_labels = einx.id("y x -> (y x)", labels)       # (R, C) → (R C,)
    sizes = np.bincount(flat_labels, minlength=count + 1)  # (R C,) → (K + 1,)
    # 1[|component k| ≥ N_min], never the background, label 0
    keep = sizes >= min_pixels                          # (K + 1,) → (K + 1,)
    keep[0] = False
    # b_n = c_n 1[|K(n)| ≥ N_min]
    return keep[labels]                                 # (K + 1,), (R, C) → (R, C)

② Strategies and ③ models

The strategies say how: LedoitWolf fits a mean and a shrunk covariance, SigmaThreshold keeps the pixels above $k$ noise levels. The model holds what was fitted: Background keeps $\gvect{\mu}$ and $\gmat{\Sigma}$, solves with $\gmat{\Sigma}$, and gives the matched filter's standard error, which the segmentation threshold uses:

• $\sigma_y \in \mathbb{R}$: the enhancement's standard error under a Gaussian background, in ppm m.
• $k \in \mathbb{R}$: the threshold, in noise levels; $\vect{c} \in \{0, 1\}^{N_\text{row} \times N_\text{col}}$: the candidate pixels, which become the plume mask $\vect{b}$ once their small components are removed.

The threshold is a one-sided test of $y_n = 0$ against $y_n > 0$. Under the Gaussian background $y_n / \sigma_y$ is standard normal where there is no methane, so each pixel is a false alarm with probability

with $\Phi$ the standard normal CDF: 0.13 % at $k = 3$, about 85 false pixels in a 256 × 256 scene, which the $N_\text{min}$ filter then removes. The candidates $\vect{c}$ become the plume mask $\vect{b}$ there.

This is where the kind axis's MatchedFilter Module splits: its state becomes the model Background, and its call becomes the algorithm matched_filter below. The split is what lets segmentation reuse $\sigma_y$.

import einx
import equinox as eqx
import jax.numpy as jnp
import lineax as lx
import numpy as np
import sklearn.covariance
from jaxtyping import Array, Bool, Float

# Shapes: W = window bands; P = valid pixels; R × C = N_row × N_col pixels


# --- ② strategies: how a model is fitted, or a rule applied ---
class LedoitWolf(eqx.Module):
    """Background mean and a shrunk covariance; swap for another estimator."""

    def __call__(self, spectra: Float[np.ndarray, "pixel band"]) -> tuple:
        """Background mean and a shrunk covariance from the valid pixels' spectra.

            μ̂ = (1 / P) ∑_n x_n
            Σ̂ = (1 − ρ) S + ρ ν I,  ρ = argmin_ρ 𝔼 ‖Σ̂ − Σ‖²_F

        with S the sample covariance and ν = tr(S) / W. Shrinking keeps Σ̂ well
        conditioned when W bands are estimated from few pixels.

        Args:
            spectra: X, one row per valid pixel.

        Returns:
            (μ̂, Σ̂).
        """
        # μ, Σ by Ledoit–Wolf shrinkage
        estimate = sklearn.covariance.LedoitWolf().fit(spectra)  # (P, W) → fitted
        return estimate.location_, estimate.covariance_  # → (W,), (W, W)


class SigmaThreshold(eqx.Module):
    """Keep pixels above k noise levels; swap for Otsu or a percentile."""

    k: float = 3.0

    def __call__(
        self, enhancement: Float[Array, "y x"], noise: float
    ) -> Bool[Array, "y x"]:
        """Pixels whose enhancement exceeds k standard errors.

            c_n = 1[y_n > k σ_y]

        A one-sided test of y_n = 0 against y_n > 0: under the Gaussian
        background a pixel is a false alarm with probability 1 − Φ(k), 0.13 %
        at k = 3, with Φ the standard normal CDF.

        Args:
            enhancement: y, ppm m.
            noise: σ_y, the matched filter's standard error, ppm m.

        Returns:
            c, the candidate pixels, before small components are removed.
        """
        # k σ_y
        level = self.k * noise                          # → ()
        # c_n = 1[y_n > k σ_y]
        return enhancement > level                      # (R, C) → (R, C)


# --- ③ models: state fitted once ---
class Background(eqx.Module):
    """The scene's background: x_n = μ + ε_n, ε_n ~ 𝒩(0, Σ), off the plume.

    `fit` estimates μ and Σ once; `solve` and `noise` use them, pure and
    jittable, for every pixel and every target t.
    """

    mean: Float[Array, " band"]
    covariance: lx.MatrixLinearOperator

    @classmethod
    def fit(
        cls, spectra: Float[np.ndarray, "pixel band"], estimator=LedoitWolf()
    ) -> "Background":
        """Estimate μ and Σ from the valid pixels' spectra, once per scene.

            (μ̂, Σ̂) = estimator(X)

        With LedoitWolf, Σ̂ = (1 − ρ) S + ρ ν I, shrunk by the ρ that minimises
        𝔼 ‖Σ̂ − Σ‖²_F. Σ̂ is kept as a PSD-tagged operator, a ⓪ structure.

        Args:
            spectra: X, one row per valid pixel.
            estimator: the strategy that gives (μ̂, Σ̂).

        Returns:
            The fitted background.
        """
        # user
        # μ, Σ from the strategy
        mean, covariance = estimator(spectra)           # (P, W) → (W,), (W, W)
        # Σ as a ⓪ structure: a PSD-tagged operator, so solves take the fast path
        operator = lx.MatrixLinearOperator(             # (W, W) → (W, W) PSD op
            jnp.asarray(covariance), lx.positive_semidefinite_tag
        )
        return cls(mean=jnp.asarray(mean), covariance=operator)

    def solve(self, vector: Float[Array, " band"]) -> Float[Array, " band"]:
        """Σ⁻¹ v, by a solve with the PSD operator, never an inverse.

        Args:
            vector: v.

        Returns:
            Σ⁻¹ v.
        """
        # Σ⁻¹ v
        solution = lx.linear_solve(self.covariance, vector)  # (W, W) op, (W,)
        return solution.value                           # → (W,)

    def noise(self, target: Float[Array, " band"]) -> Float[Array, ""]:
        """The matched filter's standard error for target t.

            σ_y² = Var(sᵀ ε) / (tᵀ s)² = sᵀ Σ s / (tᵀ s)² = 1 / (tᵀ Σ⁻¹ t)

        with s = Σ⁻¹ t and ε ~ 𝒩(0, Σ) the background's departure.

        Args:
            target: t.

        Returns:
            σ_y, in ppm m.
        """
        # Σ⁻¹ t
        direction = self.solve(target)                  # (W,) → (W,)
        # tᵀ Σ⁻¹ t
        information = einx.dot("band, band ->", target, direction)  # → ()
        # σ_y = (tᵀ Σ⁻¹ t)^(−1/2), the matched filter's standard error
        return information**-0.5                        # () → ()

④ Algorithms

Functions over models: the matched filter, the segmentation, and the emission rate from the integrated mass enhancement, with the plume length taken as the square root of the plume's area:

• $L$: the plume length in m; $U$: the wind speed at the source in m s⁻¹; $\kappa$: the mass of methane per square metre in one ppm m; $Q$: the emission rate in kg s⁻¹.

The emission rate rests on a steady-state mass balance: the plume holds what the source emitted while the wind carried it across the plume's length, so the integrated mass enhancement is

and solving for $Q$ gives the formula above. Taking $L$ as the square root of the plume's area is a convention; other lengths change $Q$ by a factor near one.

import einx
import numpy as np
from jaxtyping import Array, Bool, Float

# Shapes: W = window bands; R × C = N_row × N_col pixels


def matched_filter(
    radiance: Float[Array, "band y x"],
    background: Background,
    target: Float[Array, " band"],
) -> Float[Array, "y x"]:
    """Methane enhancement of every pixel, against a fitted background.

        y_n = argmin_y (x_n − μ − y t)ᵀ Σ⁻¹ (x_n − μ − y t) = sᵀ (x_n − μ) / (tᵀ s)

    with s = Σ⁻¹ t, from the model x_n = μ + y_n t + ε_n, ε_n ~ 𝒩(0, Σ).

    Algorithm:
        1. s = Σ⁻¹ t, through the model's solve.
        2. r_n = x_n − μ and its projection sᵀ r_n.
        3. y_n = sᵀ r_n / (tᵀ s).

    Args:
        radiance: x, the window's spectra.
        background: the fitted model, μ and Σ.
        target: t.

    Returns:
        y, in ppm m.
    """
    # user
    # s = Σ⁻¹ t, through the model
    direction = background.solve(target)                # (W,) → (W,)
    # r_n = x_n − μ
    residual = einx.subtract(                           # (W, R, C), (W,) → (W, R, C)
        "band y x, band -> band y x", radiance, background.mean
    )
    # sᵀ r_n
    projection = einx.dot("band y x, band -> y x", residual, direction)  # → (R, C)
    # tᵀ s
    normalisation = einx.dot("band, band ->", target, direction)  # → ()
    # y_n = sᵀ r_n / (tᵀ s)
    return projection / normalisation                   # (R, C) → (R, C)


def segment(
    enhancement: Float[np.ndarray, "y x"],
    noise: float,
    rule: SigmaThreshold,
    min_pixels: int,
) -> Bool[np.ndarray, "y x"]:
    """The plume mask: pixels above k σ_y, in components of at least N_min pixels.

        c_n = 1[y_n > k σ_y],  b_n = c_n 1[|K(n)| ≥ N_min]

    Algorithm:
        1. Threshold y with the rule, a strategy, so Otsu or a percentile can
           take its place.
        2. Remove the components smaller than N_min pixels, mostly noise.

    Args:
        enhancement: y, ppm m.
        noise: σ_y, ppm m.
        rule: the threshold strategy.
        min_pixels: N_min.

    Returns:
        b, the plume mask.
    """
    # user
    # c_n = 1[y_n > k σ_y], by the strategy
    above = rule(enhancement, noise)                    # (R, C) → (R, C)
    above = np.asarray(above)                           # (R, C) → (R, C)
    # user
    # b_n = c_n 1[|K(n)| ≥ N_min]
    return large_components(above, min_pixels)          # (R, C) → (R, C)


def emission_rate(
    enhancement: Float[np.ndarray, "y x"],
    plume: Bool[np.ndarray, "y x"],
    area: float,
    wind: float,
    kappa: float,
) -> float:
    """Emission rate from the integrated mass enhancement (Varon et al., 2018).

    In steady state, the methane in the plume is the source rate times the
    time the wind takes to carry it across the plume, IME = Q L / U, so

        IME = κ a ∑_n b_n y_n,  L = (a ∑_n b_n)^(1/2),  Q = U IME / L

    Algorithm:
        1. IME: the plume's enhancement summed over its pixels, in kg.
        2. L: the plume length, the square root of its area, in m.
        3. Q = U IME / L, in kg s⁻¹; zero for an empty plume.

    Args:
        enhancement: y, ppm m.
        plume: b, the plume mask.
        area: a, the pixel area, m².
        wind: U, the effective wind speed at the source, m s⁻¹.
        kappa: κ, kg m⁻² per ppm m.

    Returns:
        Q, kg s⁻¹.
    """
    # b_n y_n
    plume_enhancement = np.where(plume, enhancement, 0.0)  # (R, C) → (R, C)
    # ∑ b_n y_n, ppm m
    total = np.nansum(plume_enhancement)                # (R, C) → ()
    # IME = κ a ∑ b_n y_n, the mass of methane in the plume, kg
    ime = kappa * area * total                          # → ()
    # ∑ b_n, pixels
    plume_pixels = np.sum(plume)                        # (R, C) → ()
    # L = (a ∑ b_n)^(1/2), m
    length = np.sqrt(area * plume_pixels)               # → ()
    if length == 0:
        return 0.0
    # Q = U IME / L, kg s⁻¹
    rate = wind * ime / length                          # → ()
    return float(rate)

⑤ Steps

One Operator per stage, each a $\pi_\alpha$, a call down the ladder and an $\iota_\alpha$. Retrieval writes $\sigma_y$ into the enhancement's attributes, so segmentation reads it from its input. Quantification is two steps, because the wind lives on another grid: WindAt reads $U$ at the source's coordinates, and EmissionRate takes it as an input. In the four lifted functions $U$ was the parameter of $h[f_\text{ime}, U]$; here it is wired, as geotoolz's rule wants: a second input is an argument of _apply, never a constructor field, so a Graph can connect it.

import einx
import jax.numpy as jnp
import numpy as np
import pipekit
from geoarray import GeoArray, GeoPoints, interp_at
from geoops.carrier import wrap_like
from pipekit import Signature
from radtrax.atmosphere import air_mass_factor

# Shapes: B = 285 bands; W window bands; R × C pixels; P valid; G × H wind cells

# ch4_unit_absorption is proposed for radtrax: methane absorption per ppm m


class MethaneRetrieval(pipekit_ops.Operator):
    """Retrieval: radiance → methane enhancement, with its noise σ_y in the attrs.

        y_geo = ι_α(matched_filter(π_α(x_geo); Background.fit(X)), x_geo)
    """

    def __init__(self, *, window_nm: tuple[float, float] = (2122.0, 2488.0)) -> None:
        """Hold α, the spectral window in nm where methane absorbs strongly."""
        self.window_nm = window_nm

    def _apply(self, scene: GeoArray["band y x"]) -> GeoArray["1 y x"]:
        """Unwrap, fit the background, filter, rewrap.

            π_α(x_geo) = (x_W, u, X),  u = (M, k)
            μ̂, Σ̂ = Background.fit(X),  t = −M (k ⊙ μ̂)
            y_n = sᵀ (x_n − μ̂) / (tᵀ s),  s = Σ̂⁻¹ t,  σ_y = (tᵀ s)^(−1/2)

        Algorithm:
            1. π_α: the window's bands, u from the attrs and the wavelengths,
               and X from the valid pixels.
            2. ③ fit the model, ① build t, ④ filter.
            3. ι_α: one band in ppm m on the same grid, NaN where m = 0, and
               σ_y in 𝓜 for the steps downstream.

        Args:
            scene: x_geo, the radiance, with band wavelengths in C and the
                sun and view zenith angles in 𝓜.

        Returns:
            y_geo, the enhancement in ppm m, with noise_ppm_m = σ_y in 𝓜.
        """
        # --- π_α: the window's bands, u = (M, k) and X ---
        # today (geoarray)
        bands = scene.coord("band")
        lower, upper = self.window_nm
        names = [
            n for n, w in zip(bands.names, bands.wavelength_nm) if lower <= w <= upper
        ]
        # today (geoarray)
        swir = scene.sel(band=names)                    # (B, R, C) → (W, R, C)
        # today (radtrax)
        # M = 1 / cos ϑ_sun + 1 / cos ϑ_view, from the attrs 𝓜
        air_mass = air_mass_factor(                     # → ()
            swir.attrs["solar_zenith_deg"], swir.attrs["view_zenith_deg"]
        )
        # today (geoarray)
        wavelengths = swir.coord("band").wavelength_nm  # → W labels
        # proposed (radtrax)
        # k at the window's wavelengths
        absorption = ch4_unit_absorption(wavelengths)   # → (W,)
        # today (geoarray)
        # m: the valid pixels
        valid = swir.valid_mask()                       # (W, R, C) → (R, C)
        # X: the spectra of the valid pixels, one row each
        swir_values = np.asarray(swir)                  # (W, R, C) → (W, R, C)
        pixels = einx.id("band y x -> (y x) band", swir_values)  # → (R C, W)
        is_valid = einx.id("y x -> (y x)", valid)       # (R, C) → (R C,)
        spectra = pixels[is_valid]                      # (R C, W) → (P, W)
        # --- ③ fit the model, ① build t, ④ filter ---
        # user
        background = Background.fit(spectra)            # (P, W) → μ (W,), Σ (W, W)
        # user
        # t = −M (k ⊙ μ)
        target = target_signature(                      # (W,), (W,) → (W,)
            background.mean, jnp.asarray(absorption), air_mass
        )
        # user
        # y_n = sᵀ (x_n − μ) / (tᵀ s)
        enhancement = matched_filter(swir.data, background, target)  # → (R, C)
        # --- ι_α: one band in ppm m, NaN where m = 0, σ_y carried in 𝓜 ---
        masked = np.where(valid, enhancement, np.nan)   # (R, C) → (R, C)
        # user
        # σ_y = (tᵀ Σ⁻¹ t)^(−1/2)
        noise = background.noise(target)                # → ()
        # today (geoops)
        return wrap_like(                               # (R, C) → (1, R, C)
            swir,
            einx.id("y x -> 1 y x", masked),
            band_names=("ch4",),
            fill_value=np.nan,
            attrs={"noise_ppm_m": float(noise)},
        )

    def compute_output_signature(self, sig: Signature) -> Signature:
        """The output's dims and dtype from the input's, without data."""
        # today (pipekit)
        # (band=B, y, x) → (band=1, y, x)
        return sig.replace_dims({"band": 1})


class PlumeSegmentation(pipekit_ops.Operator):
    """Segmentation: enhancement → plume mask, at k times the retrieval's noise.

        b_geo = ι_α(segment(y, σ_y; k, N_min), y_geo)
    """

    def __init__(self, *, k: float = 3.0, min_pixels: int = 20) -> None:
        """Hold α: the threshold k in noise levels and the smallest plume N_min."""
        self.k = k
        self.min_pixels = min_pixels

    def _apply(self, enhancement: GeoArray["1 y x"]) -> GeoArray["1 y x"]:
        """Unwrap, threshold, keep the large components, rewrap.

            c_n = 1[y_n > k σ_y],  b_n = c_n 1[|K(n)| ≥ N_min]

        Algorithm:
            1. π_α: the map y and σ_y from the retrieval's attrs.
            2. ② the rule and ④ the algorithm give b.
            3. ι_α: a boolean band on the same grid.

        Args:
            enhancement: y_geo, from MethaneRetrieval.

        Returns:
            b_geo, the plume mask.
        """
        # --- π_α: the map, and σ_y from the retrieval's attrs ---
        values = einx.id("1 y x -> y x", np.asarray(enhancement))  # → (R, C)
        noise = enhancement.attrs["noise_ppm_m"]        # → ()
        # --- ② the rule, ④ the algorithm ---
        rule = SigmaThreshold(k=self.k)
        # user
        # b_n = c_n 1[|K(n)| ≥ N_min],  c_n = 1[y_n > k σ_y]
        plume = segment(values, noise, rule, self.min_pixels)  # (R, C) → (R, C)
        # --- ι_α: a boolean band on the same grid ---
        # today (geoops)
        return wrap_like(                               # (R, C) → (1, R, C)
            enhancement,
            einx.id("y x -> 1 y x", plume),
            band_names=("plume",),
            fill_value=False,
        )

    def compute_output_signature(self, sig: Signature) -> Signature:
        """The same dims as the input, as booleans."""
        # today (pipekit)
        # same grid, one boolean band
        return Signature(dims=sig.dims, dtype="bool")


class WindAt(pipekit_ops.Operator):
    """Quantification: the wind speed U at the source, read from its own grid.

        U = v(s_source), bilinear in the wind's grid
    """

    _terminal = True

    def __init__(self, *, lon: float, lat: float) -> None:
        """Hold α: the source's longitude and latitude, degrees."""
        self.lon = lon
        self.lat = lat

    def _apply(self, wind_field: GeoArray["band y x"]) -> float:
        """Read the wind at the source.

            U = ∑_corners w_c v_c,  w the bilinear weights of s_source

        Args:
            wind_field: v_geo, ERA5's 10 m wind speed on its own grid.

        Returns:
            U, m s⁻¹.
        """
        # the source, in the wind's CRS
        source_xy = np.array([[self.lon, self.lat]])    # → (1, 2)
        # today (geoarray)
        source = GeoPoints(np.zeros(1), xy=source_xy, dims="obs")  # → (1,)
        # today (geoarray)
        # U = v(s_source), bilinear
        sampled = interp_at(wind_field, source)         # (1, G, H) → (1, 1)
        values = np.asarray(sampled)                    # (1, 1) → (1, 1)
        return float(values.item())                     # (1, 1) → ()

    def compute_output_signature(self, sig: Signature) -> Signature:
        """One number, whatever the field's shape."""
        # today (pipekit)
        # one number: U, in m s⁻¹
        return Signature(dims=(), dtype="float64")


class EmissionRate(pipekit_ops.Operator):
    """Quantification: Q from the enhancement, the plume and the wind.

        Q = emission_rate(y, b; a, U, κ),  a = |det A₂ₓ₂|
    """

    _terminal = True

    def __init__(self, *, kappa: float = 6.56e-7) -> None:
        """Hold α: κ, kg m⁻² per ppm m, the near-surface value."""
        self.kappa = kappa

    def _apply(
        self, enhancement: GeoArray["1 y x"], plume: GeoArray["1 y x"], wind: float
    ) -> dict:
        """Unwrap, compute Q, return a table: the recipe's last step.

            Q = U κ a ∑_n b_n y_n / L,  L = (a ∑_n b_n)^(1/2)

        Args:
            enhancement: y_geo, ppm m.
            plume: b_geo, on the same grid.
            wind: U, m s⁻¹, from WindAt.

        Returns:
            A table, {emission_rate_kg_s: Q}; terminal, so not a carrier.
        """
        # --- π_α: u = a from A, and the two maps ---
        # user
        # a = |det A₂ₓ₂|
        area = pixel_area(enhancement.transform)        # (6,) → ()
        values = einx.id("1 y x -> y x", np.asarray(enhancement))  # → (R, C)
        mask = einx.id("1 y x -> y x", np.asarray(plume))  # → (R, C)
        # --- ④ the algorithm; terminal, so ι_α returns a table, not a carrier ---
        # user
        # Q = U κ a ∑ b_n y_n / L
        rate = emission_rate(values, mask, area, wind, self.kappa)  # → ()
        return {"emission_rate_kg_s": rate}

⑥ Recipe, and orchestration

The recipe strings the steps into one job: retrieval, segmentation, quantification. The enhancement feeds both segmentation and quantification, and the wind enters on its own grid, so the recipe is a Graph rather than a Sequential:

In symbols, four lifted functions, with the wind now an input of the last one rather than its parameter:

Orchestration runs the recipe without changing it: built in Python, summarised without data, saved as JSON, loaded back, or swept over a setting. A DVC stage, an Airflow task and a GeoApps endpoint all do one of these.

import pipekit


def plume_recipe(*, lon: float, lat: float, k: float = 3.0) -> pipekit_ops.Graph:
    """Retrieval, segmentation and quantification, wired into one Graph.

        y_geo = h[f_ret](x_geo; α₁),  b_geo = h[f_seg](y_geo; α₂)
        U = h[f_wind](v_geo; α₃),  Q = h[f_ime](y_geo, b_geo, U; α₄)

    Args:
        lon: the source's longitude, degrees, for the wind.
        lat: the source's latitude, degrees.
        k: the segmentation threshold, in noise levels.

    Returns:
        A Graph with inputs scene and wind and the output emission. It holds
        configs only, so pipekit_ops.dumps writes it and pipekit_ops.loads rebuilds it.
    """
    # today (pipekit)
    scene_in = pipekit_ops.Input("scene")
    # today (pipekit)
    wind_in = pipekit_ops.Input("wind")
    # user
    # retrieval: y_geo with σ_y
    enhancement = MethaneRetrieval()(scene_in)
    # user
    # segmentation: b_geo
    plume = PlumeSegmentation(k=k)(enhancement)
    # user
    # quantification: U, then Q
    wind = WindAt(lon=lon, lat=lat)(wind_in)
    # user
    emission = EmissionRate()(enhancement, plume, wind)
    # today (pipekit)
    return pipekit_ops.Graph(
        inputs={"scene": scene_in, "wind": wind_in}, outputs={"emission": emission}
    )


# --- orchestration: build, summarise, save, load and run; the recipe never changes ---
# user
recipe = plume_recipe(lon=-103.52, lat=31.98)
# today (pipekit)
# the signatures of the data part's carriers, with no data read
scene_sig = pipekit_ops.signature_of(scene)        # → (band=B, y=R, x=C)
# today (pipekit)
wind_sig = pipekit_ops.signature_of(wind_field)    # → (band=1, y=G, x=H)
# today (pipekit)
# a Keras-style table, from signatures alone
print(recipe.summary({"scene": scene_sig, "wind": wind_sig}))
# today (pipekit)
# the recipe as data: JSON with the steps' configs, never the arrays
recipe_json = pipekit_ops.dumps(recipe)
# today (pipekit)
# loaded back, it runs on the carriers of the data part
loaded = pipekit_ops.loads(recipe_json)
result = loaded(scene=scene, wind=wind_field)  # → {emission: table}
# user
# a sweep over the threshold, as a DVC stage or an Airflow task would run it
sweep = {
    k: plume_recipe(lon=-103.52, lat=31.98, k=k)(scene=scene, wind=wind_field)
    for k in (2.0, 3.0, 4.0)
}

recipe.summary prints the same signatures as a table, Keras-style, from signatures alone. A multi-input node shows ?: pipekit's compute_output_signature takes one input, so EmissionRate's output is not inferred.

Node                             Input sig                                                                Output sig
--------------------------------------------------------------------------------------------------------------------
[0] MethaneRetrieval(scene)      (band=285, y=256, x=256) float32                                         (band=1, y=256, x=256) float32
[1] PlumeSegmentation([0])       (band=1, y=256, x=256) float32                                           (band=1, y=256, x=256) bool
[2] WindAt(wind)                 (band=1, y=8, x=8) float64                                               () float64
[3] EmissionRate([0], [1], [2])  (band=1, y=256, x=256) float32, (band=1, y=256, x=256) bool, () float64  ?
→ emission = [3]                                                                                          ?

The recipe as data, read for people: the steps by name, each with what it takes and its config. This is the record pipekit_ops.dumps writes, with the node numbers replaced by step names; pipekit itself writes and reads the JSON below it.

inputs: [scene, wind]
steps:
  methane_retrieval:
    op: MethaneRetrieval
    in: [scene]
    window_nm: [2122.0, 2488.0]
  plume_segmentation:
    op: PlumeSegmentation
    in: [methane_retrieval]
    k: 3.0
    min_pixels: 20
  wind_at:
    op: WindAt
    in: [wind]
    lat: 31.98
    lon: -103.52
  emission_rate:
    op: EmissionRate
    in: [methane_retrieval, plume_segmentation, wind_at]
    kappa: 6.56e-07
outputs:
  emission: emission_rate

{
  "class": "Graph",
  "config": {
    "inputs": ["scene", "wind"],
    "nodes": [
      {
        "operator": 0,
        "parents": [
          { "input": "scene" }
        ]
      },
      {
        "operator": 1,
        "parents": [
          { "node": 0 }
        ]
      },
      {
        "operator": 2,
        "parents": [
          { "input": "wind" }
        ]
      },
      {
        "operator": 3,
        "parents": [
          { "node": 0 },
          { "node": 1 },
          { "node": 2 }
        ]
      }
    ],
    "operators": [
      {
        "class": "MethaneRetrieval",
        "config": {
          "window_nm": [2122.0, 2488.0]
        },
        "module": "plume"
      },
      {
        "class": "PlumeSegmentation",
        "config": {
          "k": 3.0,
          "min_pixels": 20
        },
        "module": "plume"
      },
      {
        "class": "WindAt",
        "config": {
          "lat": 31.98,
          "lon": -103.52
        },
        "module": "plume"
      },
      {
        "class": "EmissionRate",
        "config": { "kappa": 6.56e-07 },
        "module": "plume"
      }
    ],
    "outputs": {
      "emission": { "node": 3 }
    }
  },
  "module": "pipekit_ops._base.graph"
}

The six parts run as written, in order, as one module (checked on 7 Nov 2026), on a synthetic 256 × 256 scene, the size the summary shows, with a 3000 ppm m disk and a uniform 3.5 m s⁻¹ wind, with stand-ins for radtrax's air_mass_factor and the absorption lookup. The recipe returns $Q$ = 4.5 kg s⁻¹ against 4.6 kg s⁻¹ from the true plume, and the Graph loaded back from JSON returns the same. $\kappa$ = 6.56 × 10⁻⁷ kg m⁻² per ppm m is the near-surface value of geoops' ColumnToMass.