Plume operators (geo)

Plume post-processing on GeoTensors, behind the [geo] extra (pip install 'plumax[geo]'): the ten operators ported from geotoolz's plume module, as pipekit.Operators built with geoops's public carrier helpers (the pattern of radtrax's [geo] operators). Importing plumax.geo does not need the extra; the operators load on first use, and a missing extra raises an ImportError naming it.

geo

Plume operators on GeoTensors, behind the [geo] extra.

• :mod:plumax.geo.plume: the ten plume operators ported from geotoolz's plume module (#134): masks, contours, footprints, shape filters, column statistics and QND features, wind cones, the integrated mass enhancement, the cross-sectional flux and column-unit conversion.
• :mod:plumax.geo.plume_array: the array functions under them (plume_mask, wind_advection_cone, convert_column_units, ...).
• :mod:plumax.geo.vetting: PlumeVetting, :mod:plumax.vetting's candidate vetting and scene null on GeoTensors, as one Dataset (#157).

Importing this subpackage does not need the extra; the operators load on first use, and a missing extra raises an ImportError naming it.

plume — the operators

The plume operators on GeoTensors (behind the [geo] extra).

Ten operators ported from geotoolz's plume module (#134) with their behaviour unchanged; its eleventh, SBMP, moved to radtrax as radtrax.geo.MBMP. They wrap the array functions in :mod:plumax.geo.plume_array and are pipekit Operators built with geoops's public carrier helpers, as radtrax's [geo] operators are.

Per-pixel operators accept a georeader.GeoTensor or a plain np.ndarray and return the same carrier kind; geo-dependent operators (advection geometry, pixel-area maths, footprints) require a georeferenced GeoTensor and raise TypeError otherwise. The area and distance operators (PlumeFootprint, WindAdvectionCone, IMEEstimate, CrossSectionalFlux) work in metres and raise ValueError on a geographic CRS (see :func:geoops.carrier.require_projected_crs). Vector outputs are geopandas.GeoDataFrame. Algorithms follow the literature cited in :mod:plumax.geo.plume_array:

• PlumeMask / PlumeContours / PlumeFootprint: Frankenberg et al. (2016) detection and segmentation conventions.
• IMEEstimate: Varon et al. (2018) integrated mass enhancement.
• CrossSectionalFlux: Krings et al. (2011) / Varon (2018) appendix.
• PlumeShapeFilter / PlumeColumnStats / PlumeQNDFeatures: Pérez Carrasco et al. (2026), the false-plume classifier's features.

The IME and cross-section operators expect the enhancement carrier to be in kg/m^2. Use :class:ColumnToMass to convert from ppm m or mol/m^2.

Nodata: pixels that are non-finite or equal the carrier's fill_value_default (see :func:geoops.carrier.valid_pixels) never enter a threshold, statistic or sum, and are never part of a plume mask or contour. Output fill values follow the output's meaning, not the input's: boolean masks declare fill_value_default=False, float products (scores, converted columns) hold NaN at nodata pixels and declare fill_value_default=NaN, and label images use 0 (see :class:PlumeContours for its one documented exception).

PlumeMask

PlumeMask(*, threshold: ThresholdMode = 'otsu', min_area_px: int = 50, connectivity: Connectivity = 8, nbins: int = 256)

Bases: Operator

Binary plume mask from a single-band score or enhancement map.

Thresholds the input (absolute number, Otsu, or percentile) and drops connected components smaller than min_area_px pixels — the Frankenberg et al. (2016) detection convention, also used by Varon et al. (2018, 2021) for S2/AVIRIS plume identification.

Parameters:

Name	Type	Description	Default
threshold	ThresholdMode	float (absolute), "otsu", or "percentile:<p>" with p in [0, 100].	'otsu'
min_area_px	int	Minimum component size in pixels.	50
connectivity	Connectivity	4 or 8 connectivity for component labelling.	8
nbins	int	Histogram bins for the "otsu" threshold (see :func:geoops.segment.otsu_threshold); ignored by the other modes. Default 256.	256

Nodata pixels (non-finite or equal to the input's fill value) are excluded from the Otsu / percentile threshold and are always False in the mask, which declares fill_value_default=False.

Examples:

>>> mask = PlumeMask(
...     threshold="percentile:99.5", min_area_px=50,
... )(enhancement)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    threshold: ThresholdMode = "otsu",
    min_area_px: int = 50,
    connectivity: Connectivity = 8,
    nbins: int = 256,
) -> None:
    self.threshold = threshold
    self.min_area_px = min_area_px
    self.connectivity = connectivity
    self.nbins = nbins

PlumeContours

PlumeContours(*, min_area_px: int = 50, return_labels: bool = True, connectivity: Connectivity = 8)

Bases: Operator

Connected-component labelling of a plume mask.

Returns either an int32 label image (default) or a boolean mask of the retained components. Component size threshold and connectivity follow :func:label_components.

Nodata pixels (non-finite -- e.g. a NaN score -- or equal to the input's fill value) are never part of a contour. In the boolean output they are False (fill_value_default=False). In the label image they hold the input's fill value when it is a non-positive int32 value (e.g. -9999, so it cannot collide with a label); otherwise (NaN, no fill, or a positive fill) they are 0 (background), and a GeoTensor output then carries fill_value_default=0.

Examples:

>>> labels = PlumeContours(min_area_px=50)(mask)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    min_area_px: int = 50,
    return_labels: bool = True,
    connectivity: Connectivity = 8,
) -> None:
    self.min_area_px = min_area_px
    self.return_labels = return_labels
    self.connectivity = connectivity

PlumeFootprint

PlumeFootprint(*, min_area_m2: float = 500.0, simplify_tolerance: float | None = 15.0, properties: Sequence[str] | None = None, extra_properties: Sequence[Callable[..., Any]] | None = None)

Bases: Operator

Vectorize a plume mask into polygons with per-plume metadata.

The region properties come from the same machinery as :class:geoops.measure.RegionProps (defaults to :data:geoops.measure.DEFAULT_REGIONPROPS, pixel units) and the polygons from :func:georeader.vectorize.get_polygons. Returns a GeoDataFrame with one row per surviving plume: geometry, area_m2, centroid, mean_enhancement / max_enhancement (if enhancement is supplied), n_pixels, label_id, and skimage region properties such as major_axis_length, orientation, eccentricity, solidity, perimeter, bbox-*, and inertia_tensor_eigvals-*.

Nodata pixels never form a polygon (of a boolean mask, or of a label image whose fill is positive, such as a uint8 raster's 255), non-positive labels are background, and nodata enhancement pixels are ignored by the per-polygon statistics.

Geo-dependent: polygonisation needs the carrier's transform and CRS, so the input must be a georeferenced GeoTensor (plain arrays raise TypeError). Areas are in m^2 and the simplify tolerance in m, so the CRS must be projected with metre units; a geographic CRS raises ValueError (reproject with :class:geoops.geom.Reproject first).

Parameters:

Name	Type	Description	Default
min_area_m2	float	Drop polygons smaller than this area (m^2).	500.0
simplify_tolerance	float | None	Douglas-Peucker tolerance (m) for polygon simplification; None to skip.	15.0
properties	Sequence[str] | None	Region properties forwarded to regionprops_table.	None
extra_properties	Sequence[Callable[..., Any]] | None	User-defined regionprops_table callables.	None

Called as op(mask, enhancement=None). The optional enhancement carrier must sit on the mask's pixel grid (ValueError otherwise); it enables the mean_enhancement / max_enhancement statistics per polygon.

Examples:

>>> gdf = PlumeFootprint(min_area_m2=500.0)(mask, kg_m2)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    min_area_m2: float = 500.0,
    simplify_tolerance: float | None = 15.0,
    properties: Sequence[str] | None = None,
    extra_properties: Sequence[Callable[..., Any]] | None = None,
) -> None:
    self.min_area_m2 = min_area_m2
    self.simplify_tolerance = simplify_tolerance
    self.properties = tuple(
        DEFAULT_REGIONPROPS if properties is None else properties
    )
    self.extra_properties = (
        None
        if extra_properties is None
        else tuple(
            reject_config_summary(f, "extra_properties") for f in extra_properties
        )
    )

PlumeShapeFilter

PlumeShapeFilter(*, min_area_px: int = 50, min_fiber_to_major_ratio: float = 1.0, max_fiber_to_major_ratio: float = 10.0, min_fiber_width: float = 1.0, max_fiber_width: float = 100.0)

Bases: Operator

Drop instance labels that fail a fiber-shape morphology gate.

Adapted from Pérez Carrasco et al. (2026), physical_constraints. apply_physical_constraints — a post-Mask R-CNN gate that distinguishes elongated plume features from compact false positives. Per instance:

1. area = mask.sum(); drop when below min_area_px.
2. fiber_length = longest Euclidean path through the skeleton, in pixel widths (diagonal steps count sqrt(2); see :func:geoops.measure.skeleton_length).
3. major_axis = regionprops.axis_major_length — the PCA-based major axis of the equivalent ellipse. The original paper uses cv2.minAreaRect (rotated minimum-area rectangle), avoided here to keep the optional-dependency surface unchanged; the two measures agree closely for typical plume morphologies.
4. fiber_width = area / fiber_length.
5. Keep iff min_fiber_to_major_ratio <= fiber_length / major_axis <= max_fiber_to_major_ratio AND min_fiber_width <= fiber_width <= max_fiber_width.

Surviving instances keep their original integer labels; rejected instances are mapped to background (0).

Parameters:

Name	Type	Description	Default
min_area_px	int	Minimum pixel count to consider an instance.	50
min_fiber_to_major_ratio	float	Lower bound on fiber/major.	1.0
max_fiber_to_major_ratio	float	Upper bound on fiber/major.	10.0
min_fiber_width	float	Lower bound on area / fiber_length.	1.0
max_fiber_width	float	Upper bound on area / fiber_length.	100.0

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    min_area_px: int = 50,
    min_fiber_to_major_ratio: float = 1.0,
    max_fiber_to_major_ratio: float = 10.0,
    min_fiber_width: float = 1.0,
    max_fiber_width: float = 100.0,
) -> None:
    if min_area_px < 0:
        raise ValueError("min_area_px must be non-negative")
    if min_fiber_to_major_ratio < 0 or max_fiber_to_major_ratio <= 0:
        raise ValueError("fiber-to-major ratio bounds must be non-negative")
    if min_fiber_to_major_ratio > max_fiber_to_major_ratio:
        raise ValueError("min_fiber_to_major_ratio must be <= max")
    if min_fiber_width < 0 or max_fiber_width <= 0:
        raise ValueError("fiber width bounds must be non-negative")
    if min_fiber_width > max_fiber_width:
        raise ValueError("min_fiber_width must be <= max")
    self.min_area_px = int(min_area_px)
    self.min_fiber_to_major_ratio = float(min_fiber_to_major_ratio)
    self.max_fiber_to_major_ratio = float(max_fiber_to_major_ratio)
    self.min_fiber_width = float(min_fiber_width)
    self.max_fiber_width = float(max_fiber_width)

PlumeColumnStats

Bases: Operator

Per-instance in-mask vs background statistics of a column field.

Adapted from Pérez Carrasco et al. (2026), rfc_utils.xch4_metrics. For every non-zero label in the input map, returns a row with the in-mask mean, the out-of-mask (background) mean / std, and three derived metrics used by the paper's Random-Forest false-plume classifier:

• contrast = mean_plume / mean_background
• z_score = (mean_plume - mean_background) / std_background
• intensity_per_area = (sum_plume - mean_background * n_plume) / n_plume

Nodata column pixels (non-finite or equal to the column carrier's fill value) are left out of every statistic.

Called as op(labels, column): column is the column field (single-band GeoTensor or plain array) on the label map's pixel grid (ValueError otherwise).

PlumeQNDFeatures

PlumeQNDFeatures(*, degree: int = 6, eps: float = 50.0, min_samples: int = 50, perc_threshold: float = 98.0)

Bases: Operator

Per-instance Quantile-Normality-Deviation features for plume QA.

Adapted from Pérez Carrasco et al. (2026) — the feature pipeline that feeds the Random-Forest false-plume classifier (rfc_utils.get_QND_df_cms / get_plume_info). For each non-zero label in the input map this operator emits one row containing:

• qnd_column_* — polynomial coefficients of the Kolmogorov-Smirnov distance between the in-mask CDF of the column field and a Gaussian fit to the full scene's column distribution, sampled across 100 percentiles in [0.01, 0.99].
• qnd_albedo_* — same construction over the optional albedo field.
• n_clumps — DBSCAN cluster count of in-mask pixels above the perc_threshold percentile (per the paper's clump-count signal).
• contrast / z_score / intensity_per_area — the same stats that :class:PlumeColumnStats emits, included so a single call assembles the full feature vector used by the paper.

Nodata column / albedo pixels (non-finite or equal to the carrier's fill value) are left out of every feature.

Called as op(labels, column, albedo=None): column is the column-enhancement carrier and albedo an optional albedo carrier, both on the label map's pixel grid (ValueError otherwise). Without albedo the qnd_albedo_* coefficients are NaN.

Parameters:

Name	Type	Description	Default
degree	int	Degree of the QND polynomial (paper default: 6).	6
eps	float	DBSCAN eps (pixels).	50.0
min_samples	int	DBSCAN min_samples.	50
perc_threshold	float	Percentile (within the instance mask) above which pixels are eligible for clustering.	98.0

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    degree: int = 6,
    eps: float = 50.0,
    min_samples: int = 50,
    perc_threshold: float = 98.0,
) -> None:
    if degree < 1:
        raise ValueError("degree must be >= 1")
    if eps <= 0:
        raise ValueError("eps must be positive")
    if min_samples < 1:
        raise ValueError("min_samples must be >= 1")
    if not 0.0 < perc_threshold < 100.0:
        raise ValueError("perc_threshold must be in (0, 100)")
    self.degree = int(degree)
    self.eps = float(eps)
    self.min_samples = int(min_samples)
    self.perc_threshold = float(perc_threshold)

WindAdvectionCone

WindAdvectionCone(*, source: tuple[float, float], wind_u: float, wind_v: float, half_angle_deg: float = 30.0, max_distance: float = 5000.0, crs: str | None = None)

Bases: Operator

Mask likely downwind plume extent from a source and wind vector.

Rasterizes a circular sector centered at source, oriented along (wind_u, wind_v), with half-angle half_angle_deg and radius max_distance. The mask is useful as a prior for plume segmentation when the source location and wind are known (e.g. from a reanalysis product), following the practice of Varon et al. (2018, 2021).

Geo-dependent: the cone is rasterised in CRS coordinates via the carrier's transform, so the input must be a georeferenced GeoTensor (plain arrays raise TypeError). max_distance is in metres and the cone angle is measured in the map plane, so the carrier CRS must be projected with metre units; a geographic CRS raises ValueError (reproject with :class:geoops.geom.Reproject first).

Parameters:

Name	Type	Description	Default
source	tuple[float, float]	(x, y) source coordinates. If crs is supplied, source is interpreted in that CRS and reprojected to the carrier CRS (which must then be set); otherwise it is in carrier CRS coordinates (m).	required
wind_u	float	Eastward wind component (m/s).	required
wind_v	float	Northward wind component (m/s).	required
half_angle_deg	float	Half-angle of the cone in degrees, in [0, 180]; 180 is the full disc (see :func:wind_advection_cone).	30.0
max_distance	float	Radius of the cone (m).	5000.0
crs	str | None	Optional CRS of source.	None

Examples:

>>> cone = WindAdvectionCone(
...     source=(x0, y0), wind_u=3.0, wind_v=1.0,
...     half_angle_deg=30.0, max_distance=5000.0,
... )(enhancement)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    source: tuple[float, float],
    wind_u: float,
    wind_v: float,
    half_angle_deg: float = 30.0,
    max_distance: float = 5000.0,
    crs: str | None = None,
) -> None:
    self.source = as_tuple(source)
    self.wind_u = wind_u
    self.wind_v = wind_v
    self.half_angle_deg = half_angle_deg
    self.max_distance = max_distance
    self.crs = crs

IMEEstimate

IMEEstimate(*, wind_speed: float, length_method: Literal['max_axis', 'convex_hull', 'skeleton'] = 'max_axis', pixel_area_m2: float | None = None, return_uncertainty: bool = True, uncertainty_fraction: float = 0.5)

Bases: Operator

Integrated Mass Enhancement (IME) flux estimator.

Implements the IME method of Varon et al. (2018), AMT Quantifying methane point sources from fine-scale satellite observations of atmospheric methane plumes:

.. math::

\mathrm{IME} \;=\; \sum_{i \in \mathrm{plume}}
                    \Omega_i \,\Delta A_i, \qquad
Q \;=\; \frac{U_{\mathrm{eff}}\,\mathrm{IME}}{L}

where :math:\Omega_i is the column enhancement (kg/m^2) of pixel :math:i, :math:\Delta A_i is its area (m^2), :math:L is the effective plume length (m) — see :func:plume_length — and :math:U_{\mathrm{eff}} is the effective wind speed (m/s). Varon et al. (2018) calibrate :math:U_{\mathrm{eff}} against LES simulations; this operator takes the calibrated value as wind_speed.

The enhancement carrier MUST be in kg/m^2. Use :class:ColumnToMass upstream to convert from ppm m or mol/m^2: passing other units silently produces nonsense. Nodata enhancement pixels (non-finite or equal to the fill value) add nothing to the IME.

Geo-dependent: pixel area (m^2) and plume length (m) are derived from the carrier's transform, so the input must be a georeferenced GeoTensor (plain arrays raise TypeError) in a projected CRS with metre units; a geographic CRS raises ValueError (reproject with :class:geoops.geom.Reproject first).

Called as op(enhancement, plume_mask): plume_mask is a boolean carrier selecting plume pixels on the enhancement's pixel grid (same shape, transform and CRS; ValueError otherwise -- a mask on another grid would be indexed against the wrong pixels).

Parameters:

Name	Type	Description	Default
wind_speed	float	Effective wind speed :math:U_{\mathrm{eff}} in m/s.	required
length_method	Literal['max_axis', 'convex_hull', 'skeleton']	Length estimator: "max_axis", "convex_hull", or "skeleton". See :func:plume_length; "skeleton" raises ValueError for compact plumes whose skeleton collapses to a point.	'max_axis'
pixel_area_m2	float | None	Override the pixel area (m^2); default is the absolute input transform determinant in the projected CRS's metres (exact on an equal-area projection; on UTM it ignores the projection's scale factor, a < 0.1 % error near the central meridian).	None
return_uncertainty	bool	Append emission_rate_uncertainty_kg_s.	True
uncertainty_fraction	float	Fractional 1-sigma uncertainty on Q. The default 0.5 follows Varon et al. (2018) Table 3.	0.5

The call returns a dict with keys ime_kg, length_m, wind_speed_m_s, emission_rate_kg_s, and optionally emission_rate_uncertainty_kg_s.

Examples:

Estimate Q in kg/s from a ppm m enhancement map::

from plumax.geo import ColumnToMass, IMEEstimate, PlumeMask

kg_m2 = ColumnToMass(units_in="ppm_m")(enhancement)
mask = PlumeMask(threshold="percentile:99")(kg_m2)
result = IMEEstimate(
    wind_speed=3.5,
    length_method="convex_hull",
)(kg_m2, mask)
print(result["emission_rate_kg_s"])

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    wind_speed: float,
    length_method: Literal["max_axis", "convex_hull", "skeleton"] = "max_axis",
    pixel_area_m2: float | None = None,
    return_uncertainty: bool = True,
    uncertainty_fraction: float = 0.5,
) -> None:
    self.wind_speed = wind_speed
    self.length_method = length_method
    self.pixel_area_m2 = pixel_area_m2
    self.return_uncertainty = return_uncertainty
    if uncertainty_fraction < 0.0:
        raise ValueError("uncertainty_fraction must be non-negative")
    self.uncertainty_fraction = uncertainty_fraction

CrossSectionalFlux

CrossSectionalFlux(*, source: tuple[float, float], wind_u: float, wind_v: float, n_transects: int = 5, transect_spacing_m: float = 100.0)

Bases: Operator

Estimate cross-sectional flux at downwind transects.

For each transect at downwind distance :math:d_k, the operator integrates column mass across the plume:

.. math::

Q_k \;=\; U \int_{-\infty}^{\infty}
               \Omega(d_k, y) \, dy
    \;\approx\; U \sum_{i \in \mathrm{strip}_k}
               \Omega_i \,\Delta y_i

where :math:\mathrm{strip}_k is the set of plume pixels whose along-wind distance to the source falls in the half-open interval :math:[d_k - \Delta y / 2,\, d_k + \Delta y / 2), so a pixel whose centre lies on a strip boundary is counted once. This is the cross-section method of Krings et al. (2011) and Varon et al. (2018, appendix); a flat downwind profile of :math:Q_k versus :math:d_k confirms a steady-state estimate.

The enhancement carrier MUST be in kg/m^2 (see :class:ColumnToMass). Nodata enhancement pixels (non-finite or equal to the fill value) add nothing to a transect's flux. The implementation assumes (approximately) square pixels: it uses :math:\Delta y = \sqrt{\mathrm{pixel\_area}} as the across-wind pixel size. For strongly anisotropic pixel grids, reproject first.

Geo-dependent: transect placement and pixel size come from the carrier's transform and CRS, so the input must be a georeferenced GeoTensor (plain arrays raise TypeError) in a projected CRS with metre units; a geographic CRS raises ValueError (reproject with :class:geoops.geom.Reproject first).

Called as op(enhancement, plume_mask): plume_mask is a boolean plume carrier on the enhancement's pixel grid (ValueError otherwise).

Parameters:

Name	Type	Description	Default
source	tuple[float, float]	(x, y) source coordinates in the carrier CRS (m).	required
wind_u	float	Eastward wind component (m/s).	required
wind_v	float	Northward wind component (m/s).	required
n_transects	int	Number of downwind transects to evaluate.	5
transect_spacing_m	float	Along-wind distance between transects (m).	100.0

The call returns a geopandas.GeoDataFrame with one row per transect: geometry (across-wind LineString), transect_id, distance_m, flux_kg_s, n_pixels.

Examples:

Evaluate fluxes at 100 m, 200 m, 300 m downwind::

gdf = CrossSectionalFlux(
    source=(x0, y0),
    wind_u=3.0, wind_v=1.0,
    n_transects=3, transect_spacing_m=100.0,
)(kg_m2, mask)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    source: tuple[float, float],
    wind_u: float,
    wind_v: float,
    n_transects: int = 5,
    transect_spacing_m: float = 100.0,
) -> None:
    self.source = as_tuple(source)
    self.wind_u = wind_u
    self.wind_v = wind_v
    self.n_transects = n_transects
    self.transect_spacing_m = transect_spacing_m

ColumnToMass

ColumnToMass(*, gas: str = 'CH4', units_in: ColumnUnit = 'ppm_m', units_out: ColumnUnit = 'kg_m2')

Bases: Operator

Convert column enhancement units (ppm m / mol m^-2 / kg m^-2).

Thin wrapper over :func:convert_column_units. The ppm m conversion assumes a standard molar volume (298.15 K, 1 atm); see that function's docstring for the exact relation and caveats. Nodata pixels (non-finite or equal to the fill value) are not converted: they hold NaN and the output declares fill_value_default=NaN (a converted column is a new quantity, so the input's fill -- e.g. -9999 ppm m -- is not carried over).

Parameters:

Name	Type	Description	Default
gas	str	"CH4" (default) or "CO2".	'CH4'
units_in	ColumnUnit	Units of the input carrier.	'ppm_m'
units_out	ColumnUnit	Desired output units.	'kg_m2'

Examples:

>>> kg_m2 = ColumnToMass(
...     gas="CH4", units_in="ppm_m", units_out="kg_m2",
... )(enhancement)

Source code in src/plumax/geo/plume.py

def __init__(
    self,
    *,
    gas: str = "CH4",
    units_in: ColumnUnit = "ppm_m",
    units_out: ColumnUnit = "kg_m2",
) -> None:
    self.gas = gas
    self.units_in = units_in
    self.units_out = units_out

vetting — plume vetting on carriers

PlumeVetting runs plume vetting (vet_candidate and the scene null) on every candidate of a scene and returns the scores as one xarray.Dataset; the decision rule is left to the caller.

Plume vetting on GeoTensors (behind the [geo] extra): one Dataset per scene.

PlumeVetting is the L2 operator over the array-level vetting in :mod:plumax.vetting (#157): it reads a radiance GeoTensor, its matched-filter GeoTensor and a stack of candidate masks, runs vet_candidate on each candidate and the scene null (null_masks, null_pvalue) beside it, and returns the scores as one xarray.Dataset. It returns scores, not a verdict: the decision rule is #158's.

PlumeVetting

PlumeVetting(window: AbsorptionWindow, table: Float[ArrayLike, 'D Bf'], delta_grid: Float[ArrayLike, ' D'], *, n_null: int = 100, key: Array | None = None, k: int = 40, keep: float = 0.5, degree: int = 2, box: int = 200)

Bases: Operator

Vet every candidate of a scene and its scene null; scores as a Dataset.

Called as PlumeVetting(window, table, delta_grid)(cube, mf, masks). For each candidate M_c of masks, with every candidate kept out of the background pool,

(ΔX̂_c, D_norm,c) = vet_candidate(L, MF, M_c), and

p̂_c = null_pvalue(D_norm,c, {D_norm(M_s)}) over n_null copies M_s of M_c from null_masks, each vetted as M_c is; a copy that cannot be scored (radtrax's DegenerateInputError) scores NaN and null_pvalue leaves it out.

The cube is (B, H, W) with its band centres in nm in attrs["wavelengths"]; mf is (H, W) or (1, H, W) on the cube's pixel grid, in mol/m² (the units of window.mf_background_mol_m2). Nodata pixels of either carrier (non-finite, or equal to its fill_value_default) are dropped from every candidate and kept out of the pool and the null copies. The candidates' random draws split key in candidate order, so the same key gives the same p_null.

The Dataset has dims candidate (C) and band (the window's Bf fitting bands, coordinate wavelength_nm) and holds delta_x_mol_m2 (ΔX̂), d_norm, p_null (NaN when n_null = 0) and n_pairs per candidate; ratio (R), continuum (P(λ; ĉ)) and fit (P(λ; ĉ) E(λ; ΔX̂), the fitted ratio) per candidate and band; and attributes window (the window's fields as JSON) and key (the key's raw data).

Parameters:

Name	Type	Description	Default
window	AbsorptionWindow	The gas's AbsorptionWindow (radtrax.retrieval.absorption).	required
table	Float[ArrayLike, 'D Bf']	E(λ; ΔX_d) on the window's fitting bands, (D, Bf), from radtrax's transmittance_table.	required
delta_grid	Float[ArrayLike, ' D']	The columns ΔX_d of the table's rows [mol/m²], (D,), starting at 0.	required
n_null	int	Number of null copies per candidate; 0 skips the null.	100
key	Array | None	A jax.random key for the null copies; None uses jax.random.key(0).	None
k	int	Number of seeds before growth (vet_candidate).	40
keep	float	Fraction of the pairs kept, in (0, 1].	0.5
degree	int	The continuum's Legendre degree.	2
box	int	Margin, in pixels, around the seeds' bounding box that the pool is drawn from.	200

Raises:

Type	Description
ParameterError	If n_null is negative.
ShapeError	If the cube is not (B, H, W), has no attrs["wavelengths"] with B entries, table does not have one column per fitting band, or masks is not (C, H, W) on the cube's grid; and whatever vet_candidate raises on a candidate.

Examples:

A pipeline sketch, after radtrax's matched filter::

from plumax.geo import PlumeVetting

scores = PlumeVetting(window, table, delta_grid, n_null=100)(
    cube, mf, masks
)
scores["d_norm"], scores["p_null"]

Source code in src/plumax/geo/vetting.py

def __init__(
    self,
    window: AbsorptionWindow,
    table: Float[ArrayLike, "D Bf"],
    delta_grid: Float[ArrayLike, " D"],
    *,
    n_null: int = 100,
    key: Array | None = None,
    k: int = 40,
    keep: float = 0.5,
    degree: int = 2,
    box: int = 200,
) -> None:
    if n_null < 0:
        raise ParameterError(f"PlumeVetting: n_null must be >= 0, got {n_null}")
    self.window = window
    self.table = table
    self.delta_grid = delta_grid
    self.n_null = int(n_null)
    self.key = jax.random.key(0) if key is None else key
    self.k = k
    self.keep = keep
    self.degree = degree
    self.box = box

get_config

get_config() -> dict[str, Any]

The JSON-safe configuration, without the arrays (table, grid, key).

Source code in src/plumax/geo/vetting.py

def get_config(self) -> dict[str, Any]:
    """The JSON-safe configuration, without the arrays (table, grid, key)."""
    return {
        "window": json.loads(json.dumps(dataclasses.asdict(self.window))),
        "n_null": self.n_null,
        "k": self.k,
        "keep": self.keep,
        "degree": self.degree,
        "box": self.box,
    }

from_config classmethod

from_config(config: dict[str, Any], *, table: Float[ArrayLike, 'D Bf'], delta_grid: Float[ArrayLike, ' D'], key: Array | None = None) -> PlumeVetting

Rebuild from get_config's output and the arrays it leaves out.

Source code in src/plumax/geo/vetting.py

@classmethod
def from_config(
    cls,
    config: dict[str, Any],
    *,
    table: Float[ArrayLike, "D Bf"],
    delta_grid: Float[ArrayLike, " D"],
    key: Array | None = None,
) -> PlumeVetting:
    """Rebuild from `get_config`'s output and the arrays it leaves out."""
    config = dict(config)
    window = dict(config.pop("window"))
    for name in ("fit_nm", "absorbing_nm"):
        window[name] = tuple(tuple(r) for r in window[name])
    return cls(
        AbsorptionWindow(**window),
        table,
        delta_grid,
        key=key,
        **config,
    )

plume_array — the array functions

Plume masks, geometry and column units on NumPy arrays (behind [geo]).

The array functions under the operators in :mod:plumax.geo.plume, ported from geotoolz's plume module (#134) with their behaviour unchanged. They take plain arrays and an affine transform, no carrier, but build on geoops's public mask, measure and threshold helpers and on shapely, so they need the [geo] extra too. Algorithms and unit conventions follow the trace-gas plume literature cited per function:

• Varon et al. (2018), AMT: the integrated mass enhancement (IME) method.
• Frankenberg et al. (2016), PNAS: airborne CH4 plume detection.
• Thompson et al. (2015), AMT, and Foote et al. (2020), TGRS: the matched-filter CH4 retrievals whose columns these functions convert.

plume_mask

plume_mask(values: Num[ndarray, 'h w'] | Num[ndarray, '1 h w'], *, threshold: ThresholdMode = 'otsu', min_area_px: int = 50, connectivity: Connectivity = 8, nbins: int = 256) -> Bool[np.ndarray, 'h w']

Threshold an enhancement map and remove small connected components.

Parameters:

Name	Type	Description	Default
values	Num[ndarray, 'h w'] | Num[ndarray, '1 h w']	Single-band enhancement map, (H, W) or (1, H, W).	required
threshold	ThresholdMode	Absolute number, "otsu", or "percentile:<p>"; see :func:geoops.segment.resolve_threshold.	'otsu'
min_area_px	int	Minimum connected-component size in pixels.	50
connectivity	Connectivity	4 or 8 connectivity for component labelling.	8
nbins	int	Histogram bins for the "otsu" threshold. Default 256.	256

Returns:

Type	Description
Bool[ndarray, 'h w']	Boolean (H, W) mask of the surviving plume pixels.

Source code in src/plumax/geo/plume_array.py

def plume_mask(
    values: Num[np.ndarray, "h w"] | Num[np.ndarray, "1 h w"],
    *,
    threshold: ThresholdMode = "otsu",
    min_area_px: int = 50,
    connectivity: Connectivity = 8,
    nbins: int = 256,
) -> Bool[np.ndarray, "h w"]:
    """Threshold an enhancement map and remove small connected components.

    Args:
        values: Single-band enhancement map, ``(H, W)`` or ``(1, H, W)``.
        threshold: Absolute number, ``"otsu"``, or ``"percentile:<p>"``;
            see :func:`geoops.segment.resolve_threshold`.
        min_area_px: Minimum connected-component size in pixels.
        connectivity: 4 or 8 connectivity for component labelling.
        nbins: Histogram bins for the ``"otsu"`` threshold. Default ``256``.

    Returns:
        Boolean ``(H, W)`` mask of the surviving plume pixels.
    """
    arr = squeeze_single_band(values)
    raw = threshold_mask(arr, threshold, nbins=nbins)
    return label_components(raw, min_area_px=min_area_px, connectivity=connectivity) > 0

wind_advection_cone

wind_advection_cone(shape: tuple[int, int], transform: Any, *, source: tuple[float, float], wind_u: float, wind_v: float, half_angle_deg: float = 30.0, max_distance: float = 5000.0) -> Bool[np.ndarray, 'h w']

Rasterize an analytical downwind sector mask.

Marks pixels whose center lies within max_distance of source and whose bearing from the source is within half_angle_deg of the wind direction (wind_u, wind_v) — a geometric prior for where an advected plume can be. The source pixel itself is included.

With d = (x − x₀, y − y₀) the offset of a pixel center from the source and ŵ = (u, v) / ‖(u, v)‖ the unit wind vector, a pixel is inside when

‖d‖ ≤ max_distance  and  d·ŵ ≥ ‖d‖ · cos(half_angle_deg)

i.e. the angle between d and ŵ is at most half_angle_deg. For half_angle_deg ≤ 90 the angle test alone keeps the sector in the downwind half-plane (cos ≥ 0 ⇒ d·ŵ ≥ 0); half_angle_deg > 90 widens it past the crosswind line, and half_angle_deg = 180 is the full disc of radius max_distance, upwind pixels included.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Raster shape (H, W).	required
transform	Any	Affine-like geotransform of the raster; its CRS units must match source and max_distance.	required
source	tuple[float, float]	(x, y) source coordinates in CRS units.	required
wind_u	float	Eastward wind component.	required
wind_v	float	Northward wind component.	required
half_angle_deg	float	Half-angle of the sector in degrees, in [0, 180].	30.0
max_distance	float	Sector radius in CRS units.	5000.0

Returns:

Type	Description
Bool[ndarray, 'h w']	Boolean (H, W) mask, True inside the downwind sector.

Raises:

Type	Description
ValueError	If the wind vector is zero, half_angle_deg is outside [0, 180], or max_distance is not positive.

Source code in src/plumax/geo/plume_array.py

def wind_advection_cone(
    shape: tuple[int, int],
    transform: Any,
    *,
    source: tuple[float, float],
    wind_u: float,
    wind_v: float,
    half_angle_deg: float = 30.0,
    max_distance: float = 5000.0,
) -> Bool[np.ndarray, "h w"]:
    """Rasterize an analytical downwind sector mask.

    Marks pixels whose center lies within ``max_distance`` of ``source``
    and whose bearing from the source is within ``half_angle_deg`` of the
    wind direction ``(wind_u, wind_v)`` — a geometric prior for where an
    advected plume can be. The source pixel itself is included.

    With ``d = (x − x₀, y − y₀)`` the offset of a pixel center from the
    source and ``ŵ = (u, v) / ‖(u, v)‖`` the unit wind vector, a pixel is
    inside when

        ‖d‖ ≤ max_distance  and  d·ŵ ≥ ‖d‖ · cos(half_angle_deg)

    i.e. the angle between ``d`` and ``ŵ`` is at most ``half_angle_deg``.
    For ``half_angle_deg ≤ 90`` the angle test alone keeps the sector in
    the downwind half-plane (``cos ≥ 0 ⇒ d·ŵ ≥ 0``); ``half_angle_deg > 90``
    widens it past the crosswind line, and ``half_angle_deg = 180`` is the
    full disc of radius ``max_distance``, upwind pixels included.

    Args:
        shape: Raster shape ``(H, W)``.
        transform: Affine-like geotransform of the raster; its CRS units
            must match ``source`` and ``max_distance``.
        source: ``(x, y)`` source coordinates in CRS units.
        wind_u: Eastward wind component.
        wind_v: Northward wind component.
        half_angle_deg: Half-angle of the sector in degrees, in [0, 180].
        max_distance: Sector radius in CRS units.

    Returns:
        Boolean ``(H, W)`` mask, True inside the downwind sector.

    Raises:
        ValueError: If the wind vector is zero, ``half_angle_deg`` is
            outside [0, 180], or ``max_distance`` is not positive.
    """
    wind_norm = float(np.hypot(wind_u, wind_v))
    if wind_norm == 0.0:
        raise ValueError("wind vector must be non-zero")
    if not 0.0 <= half_angle_deg <= 180.0:
        raise ValueError("half_angle_deg must be in [0, 180]")
    if max_distance <= 0.0:
        raise ValueError("max_distance must be positive")

    xs, ys = pixel_centers(shape, transform)
    dx = xs - source[0]
    dy = ys - source[1]
    distances = np.hypot(dx, dy)
    projection = (dx * wind_u + dy * wind_v) / wind_norm
    cos_angle = np.divide(
        projection,
        distances,
        out=np.ones_like(distances, dtype=float),
        where=distances > 0,
    )
    # Clip rounding overshoot so half_angle_deg=180 (cos = −1) is the full disc.
    cos_angle = np.clip(cos_angle, -1.0, 1.0)
    min_cos = np.cos(np.deg2rad(half_angle_deg))
    return (distances <= max_distance) & (cos_angle >= min_cos)

convert_column_units

convert_column_units(values: Num[ndarray, '*dims'], *, gas: str = 'CH4', units_in: ColumnUnit = 'ppm_m', units_out: ColumnUnit = 'kg_m2') -> Float[np.ndarray, '*dims']

Convert column enhancement among ppm m, mol/m^2, and kg/m^2.

The conversions are

.. math::

\Omega_{\mathrm{mol/m^2}} \;=\; \frac{X \cdot 10^{-6}}{V_m}, \qquad
\Omega_{\mathrm{kg/m^2}}  \;=\; M_{\mathrm{gas}}
                                \cdot \Omega_{\mathrm{mol/m^2}}

where X is the column-integrated volume mixing ratio in ppm m, :math:V_m = 0.024465 m^3/mol is the standard molar volume of an ideal gas at 298.15 K and 1 atm, and :math:M_{\mathrm{gas}} is the molar mass (CH4: 0.01604 kg/mol; CO2: 0.04401 kg/mol).

The ppm m branch therefore assumes near-surface conditions. Pass mol_m2 or kg_m2 inputs when retrieval-specific pressure and temperature corrections have already been applied upstream (e.g. by a matched-filter retrieval per Thompson 2015 / Foote 2020).

Source code in src/plumax/geo/plume_array.py

def convert_column_units(
    values: Num[np.ndarray, "*dims"],
    *,
    gas: str = "CH4",
    units_in: ColumnUnit = "ppm_m",
    units_out: ColumnUnit = "kg_m2",
) -> Float[np.ndarray, "*dims"]:
    r"""Convert column enhancement among ppm m, mol/m^2, and kg/m^2.

    The conversions are

    .. math::

        \Omega_{\mathrm{mol/m^2}} \;=\; \frac{X \cdot 10^{-6}}{V_m}, \qquad
        \Omega_{\mathrm{kg/m^2}}  \;=\; M_{\mathrm{gas}}
                                        \cdot \Omega_{\mathrm{mol/m^2}}

    where ``X`` is the column-integrated volume mixing ratio in ppm m,
    :math:`V_m = 0.024465` m^3/mol is the standard molar volume of an
    ideal gas at 298.15 K and 1 atm, and :math:`M_{\mathrm{gas}}` is the
    molar mass (CH4: 0.01604 kg/mol; CO2: 0.04401 kg/mol).

    The ppm m branch therefore assumes near-surface conditions. Pass
    ``mol_m2`` or ``kg_m2`` inputs when retrieval-specific pressure and
    temperature corrections have already been applied upstream (e.g. by
    a matched-filter retrieval per Thompson 2015 / Foote 2020).
    """
    gas_key = gas.upper()
    if gas_key not in MOLAR_MASS_KG_PER_MOL:
        supported = ", ".join(sorted(MOLAR_MASS_KG_PER_MOL))
        raise ValueError(f"unsupported gas {gas!r}; expected one of {supported}")
    molar_mass = MOLAR_MASS_KG_PER_MOL[gas_key]

    arr = np.asarray(values, dtype=float)
    if units_in == "ppm_m":
        # 1e-6 converts ppm to a fraction; molar volume assumes standard air.
        mol_m2 = arr * 1e-6 / STANDARD_MOLAR_VOLUME_M3_PER_MOL
    elif units_in == "mol_m2":
        mol_m2 = arr
    elif units_in == "kg_m2":
        mol_m2 = arr / molar_mass
    else:
        raise ValueError("units_in must be 'ppm_m', 'mol_m2', or 'kg_m2'")

    if units_out == "ppm_m":
        return mol_m2 * STANDARD_MOLAR_VOLUME_M3_PER_MOL * 1e6
    if units_out == "mol_m2":
        return mol_m2
    if units_out == "kg_m2":
        return mol_m2 * molar_mass
    raise ValueError("units_out must be 'ppm_m', 'mol_m2', or 'kg_m2'")

plume_length

plume_length(mask: Bool[ndarray, 'h w'], transform: Any, *, method: Literal['max_axis', 'convex_hull', 'skeleton'] = 'max_axis') -> float

Estimate plume length L from active pixel centers.

L is the effective length used by Varon et al. (2018) in the IME method Q = U_eff * IME / L. Three estimators are supported:

• "max_axis": maximum pairwise distance between active pixels, taken over the vertices of their convex hull (the same value as over every pair, in memory that does not grow with N²). Fast and robust; matches the original Varon 2018 definition for reasonably linear plumes.
• "convex_hull": max diameter of the convex hull of the active pixels. For degenerate hulls (single point, collinear points) we fall back to the point-set diameter so the result is always well defined.
• "skeleton": length of the plume centreline, via :func:geoops.measure.skeleton_length with the transform as the pixel step. The mask is skeletonised and L is the longest geodesic path through the 8-connected skeleton, with each step weighted by its Euclidean length in CRS units (so diagonal steps and anisotropic pixels are measured correctly). Better than the chord for curved plumes. The path runs between skeleton pixel centres, so for thick plumes it is shorter than the mask extent by roughly the plume width.

Lengths are in CRS units (m for a projected metric CRS) and measured between pixel centres.

A single-pixel mask returns one linear pixel size, sqrt(pixel_area), for every method: it has no geometry to measure, and the one-pixel proxy is the documented convention.

Parameters:

Name	Type	Description	Default
mask	Bool[ndarray, 'h w']	Boolean (H, W) plume mask.	required
transform	Any	Affine-like geotransform (a, b, d, e coefficients) mapping pixel indices to CRS coordinates.	required
method	Literal['max_axis', 'convex_hull', 'skeleton']	Length estimator (see above).	'max_axis'

Returns:

Type	Description
float	The plume length in CRS units; 0.0 for an empty mask.

Raises:

Type	Description
ValueError	If method is unknown, or if method="skeleton" and a multi-pixel mask yields L == 0 (every skeleton component collapses to a single pixel, e.g. a compact blob a few pixels across or scattered isolated pixels). A one-pixel proxy would badly under-estimate L (and over-estimate Q) there, so use "max_axis" or "convex_hull".

Source code in src/plumax/geo/plume_array.py

def plume_length(
    mask: Bool[np.ndarray, "h w"],
    transform: Any,
    *,
    method: Literal["max_axis", "convex_hull", "skeleton"] = "max_axis",
) -> float:
    """Estimate plume length ``L`` from active pixel centers.

    ``L`` is the effective length used by Varon et al. (2018) in the IME
    method ``Q = U_eff * IME / L``. Three estimators are supported:

    - ``"max_axis"``: maximum pairwise distance between active pixels,
      taken over the vertices of their convex hull (the same value as over
      every pair, in memory that does not grow with ``N²``). Fast and robust; matches the original Varon 2018 definition for
      reasonably linear plumes.
    - ``"convex_hull"``: max diameter of the convex hull of the active
      pixels. For degenerate hulls (single point, collinear points) we
      fall back to the point-set diameter so the result is always well
      defined.
    - ``"skeleton"``: length of the plume centreline, via
      :func:`geoops.measure.skeleton_length` with the transform as the
      pixel step. The mask is skeletonised and ``L`` is the longest
      geodesic path through the 8-connected skeleton,
      with each step weighted by its Euclidean length in CRS units
      (so diagonal steps and anisotropic pixels are measured
      correctly). Better than the chord for curved plumes. The path
      runs between skeleton pixel centres, so for thick plumes it is
      shorter than the mask extent by roughly the plume width.

    Lengths are in CRS units (m for a projected metric CRS) and measured
    between pixel centres.

    A single-pixel mask returns one linear pixel size,
    ``sqrt(pixel_area)``, for every method: it has no geometry to
    measure, and the one-pixel proxy is the documented convention.

    Args:
        mask: Boolean ``(H, W)`` plume mask.
        transform: Affine-like geotransform (``a``, ``b``, ``d``, ``e``
            coefficients) mapping pixel indices to CRS coordinates.
        method: Length estimator (see above).

    Returns:
        The plume length in CRS units; ``0.0`` for an empty mask.

    Raises:
        ValueError: If ``method`` is unknown, or if ``method="skeleton"``
            and a multi-pixel mask yields ``L == 0`` (every skeleton
            component collapses to a single pixel, e.g. a compact blob a
            few pixels across or scattered isolated pixels). A one-pixel
            proxy would badly under-estimate ``L`` (and over-estimate
            ``Q``) there, so use ``"max_axis"`` or ``"convex_hull"``.
    """
    active = np.asarray(mask, dtype=bool)
    if not active.any():
        return 0.0
    xs, ys = pixel_centers(active.shape, transform)
    points = np.column_stack([xs[active], ys[active]])
    if points.shape[0] == 1:
        # Single pixel: use one linear pixel size as the length proxy.
        return float(np.sqrt(pixel_area(transform)))
    if method == "max_axis":
        return _diameter(points)
    if method == "convex_hull":
        hull = MultiPoint(points).convex_hull
        # For polygonal hulls use the exterior ring; for degenerate
        # hulls (Point, LineString) shapely exposes ``.coords`` instead.
        exterior = getattr(hull, "exterior", None)
        if exterior is not None:
            coords = np.asarray(exterior.coords)
        elif hasattr(hull, "coords"):
            coords = np.asarray(hull.coords)
        else:
            coords = points
        diff = coords[:, None, :] - coords[None, :, :]
        return float(np.sqrt(np.max(np.sum(diff**2, axis=-1))))
    if method == "skeleton":
        length = skeleton_length(active, step=transform)
        if length <= 0.0:
            raise ValueError(
                f"skeleton plume length is 0 for a non-empty {int(active.sum())}-"
                "pixel mask: every skeleton component collapses to a single "
                "pixel (compact blob or isolated pixels). Use "
                "length_method='max_axis' or 'convex_hull' for such plumes."
            )
        return length
    raise ValueError("length_method must be 'max_axis', 'convex_hull', or 'skeleton'")

pixel_area

pixel_area(transform: Any) -> float

Return pixel area from an affine-like transform determinant.

Parameters:

Name	Type	Description	Default
transform	Any	Affine-like object exposing a, b, d, e coefficients (e.g. rasterio.Affine).	required

Returns:

Type	Description
float	|a*e - b*d| — the pixel area in squared CRS units (m^2 for a
float	projected metric CRS).

Source code in src/plumax/geo/plume_array.py

def pixel_area(transform: Any) -> float:
    """Return pixel area from an affine-like transform determinant.

    Args:
        transform: Affine-like object exposing ``a``, ``b``, ``d``, ``e``
            coefficients (e.g. ``rasterio.Affine``).

    Returns:
        ``|a*e - b*d|`` — the pixel area in squared CRS units (m^2 for a
        projected metric CRS).
    """
    return float(abs(transform.a * transform.e - transform.b * transform.d))

pixel_centers

pixel_centers(shape: tuple[int, int], transform: Any) -> tuple[Float[np.ndarray, 'h w'], Float[np.ndarray, 'h w']]

Return x/y coordinate grids for pixel centers.

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Raster shape (H, W).	required
transform	Any	Affine-like geotransform mapping (col, row) pixel indices to CRS coordinates.	required

Returns:

Type	Description
Float[ndarray, 'h w']	Tuple (xs, ys) of (H, W) float arrays holding the CRS
Float[ndarray, 'h w']	coordinates of each pixel center.

Source code in src/plumax/geo/plume_array.py

def pixel_centers(
    shape: tuple[int, int], transform: Any
) -> tuple[Float[np.ndarray, "h w"], Float[np.ndarray, "h w"]]:
    """Return x/y coordinate grids for pixel centers.

    Args:
        shape: Raster shape ``(H, W)``.
        transform: Affine-like geotransform mapping ``(col, row)`` pixel
            indices to CRS coordinates.

    Returns:
        Tuple ``(xs, ys)`` of ``(H, W)`` float arrays holding the CRS
        coordinates of each pixel center.
    """
    rows, cols = np.indices(shape, dtype=float)
    return pixel_xy(transform, rows, cols)

squeeze_single_band

squeeze_single_band(values: Shaped[ndarray, 'h w'] | Shaped[ndarray, '1 h w']) -> Shaped[np.ndarray, 'h w']

Return a 2-D plume map from a 2-D or singleton-band array.

Thin delegate to :func:geoops.carrier.single_band, kept under geotoolz's historical public name.

Parameters:

Name	Type	Description	Default
values	Shaped[ndarray, 'h w'] | Shaped[ndarray, '1 h w']	A (H, W) array or a (1, H, W) singleton-band cube. Any array-like is accepted.	required

Returns:

Type	Description
Shaped[ndarray, 'h w']	The (H, W) array. No copy is made for ndarray input.

Raises:

Type	Description
ValueError	If values is neither (H, W) nor (1, H, W).

Source code in src/plumax/geo/plume_array.py

def squeeze_single_band(
    values: Shaped[np.ndarray, "h w"] | Shaped[np.ndarray, "1 h w"],
) -> Shaped[np.ndarray, "h w"]:
    """Return a 2-D plume map from a 2-D or singleton-band array.

    Thin delegate to :func:`geoops.carrier.single_band`, kept under
    geotoolz's historical public name.

    Args:
        values: A ``(H, W)`` array or a ``(1, H, W)`` singleton-band
            cube. Any array-like is accepted.

    Returns:
        The ``(H, W)`` array. No copy is made for ndarray input.

    Raises:
        ValueError: If ``values`` is neither ``(H, W)`` nor ``(1, H, W)``.
    """
    return single_band(values, name="plume")