Quantifying Methane Plume Emission Rates in Python #

Detecting a methane plume is the easy half. Turning that plume into a number an auditor will accept — kilograms per hour, with an interval you can defend — is where most implementations quietly go wrong. This guide implements the two mainstream quantification methods end to end, and sits within methane plume detection from hyperspectral imagery, part of the satellite imagery processing stack. It assumes you already have a per-pixel enhancement field in ppm·m and a plume mask; if you do not, start with the retrieval architecture in the parent topic.

The two methods answer the same question differently. Integrated mass enhancement (IME) measures how much excess methane is sitting in the plume right now and divides by how long it takes the wind to clear it. Cross-sectional flux (CSF) measures how much methane crosses a line drawn perpendicular to the plume axis per unit time. IME is robust for compact plumes and tolerant of a ragged mask edge; CSF is better for long, well-developed plumes and gives you several independent estimates from one image. Production pipelines compute both and reconcile them, because agreement between two structurally different estimators is much stronger evidence than either alone.

Root Cause Analysis #

Flux is not a measured quantity. What the instrument measures is a column enhancement — excess methane along the light path, in ppm·m — over a set of pixels at one instant. Converting that instantaneous, static picture into a rate requires a transport model, and every quantification error traces back to an assumption inside that model.

The first assumption is steady state: that the plume you see was produced by a source emitting at a constant rate for at least as long as the plume’s transit time. For a plume 500 m long in a 3 m s⁻¹ wind, that is roughly three minutes — usually fine. For a kilometre-scale plume in light wind it can be fifteen minutes or more, during which a compressor blowdown may have started and stopped. When steady state is violated, IME reports a time-average weighted by an unknown history, and the honest response is to flag the estimate rather than to report it as an instantaneous rate.

The second assumption is the effective wind speed. The wind that transports the plume is not the wind reported at 10 m by a reanalysis product; it is a vertically integrated speed over the plume’s actual depth, in a boundary layer whose shear depends on stability and roughness. The literature standard is an empirical scaling from the 10 m wind, of the form Uₑ = a·U₁₀ + b, fitted against large-eddy simulations of plumes with known source rates. The coefficients differ by instrument and pixel size because the fit absorbs the retrieval’s own smoothing. Borrowing a scaling fitted for one instrument and applying it to another is a real and common error, worth 10–20% bias.

The third assumption is the plume mask boundary. IME sums enhancement over the mask, so any threshold choice trades a low bias against including noise. Set the threshold too high and you truncate the plume’s diffuse tail, losing mass; too low and you sum noise that scales with the mask’s area. The practical resolution is not to find the perfect threshold but to compute IME over a range of thresholds and report the plateau — a stable region where the estimate is insensitive to the choice. Absence of a plateau is itself diagnostic: it means the plume is not separable from the background.

Diagnostic Pipeline / Pre-Flight Validation #

Before quantifying, check that the plume is quantifiable. The pre-flight below rejects three conditions that make any flux estimate meaningless: an unresolved plume smaller than a few pixels, a mask touching the scene edge (so an unknown fraction is missing), and a plume whose principal axis disagrees with the wind bearing by more than a tolerance, which usually means the mask has merged two sources or captured an artefact.

import numpy as np
import structlog
from scipy import ndimage

log = structlog.get_logger()

MIN_PIXELS = 9
MAX_AXIS_WIND_DISAGREEMENT_DEG = 45.0


def preflight(mask: np.ndarray, enhancement: np.ndarray, wind_bearing_deg: float,
              pixel_size_m: float) -> dict:
    """Decide whether this plume can be quantified at all. Returns a verdict dict;
    a False verdict is a legitimate output, not an error."""
    pixels = int(mask.sum())
    if pixels < MIN_PIXELS:
        return {"quantifiable": False, "reason": "unresolved", "pixels": pixels}

    ys, xs = np.nonzero(mask)
    if ys.min() == 0 or xs.min() == 0 or ys.max() == mask.shape[0] - 1 \
            or xs.max() == mask.shape[1] - 1:
        # Truncated plumes give a guaranteed low bias of unknown size.
        return {"quantifiable": False, "reason": "touches_scene_edge", "pixels": pixels}

    # Principal axis via the mass-weighted second moment of the enhancement.
    weights = enhancement[ys, xs]
    cy, cx = np.average(ys, weights=weights), np.average(xs, weights=weights)
    cov = np.cov(np.c_[xs - cx, ys - cy].T, aweights=weights)
    eigvals, eigvecs = np.linalg.eigh(cov)
    axis = eigvecs[:, np.argmax(eigvals)]
    axis_bearing = (np.degrees(np.arctan2(axis[0], -axis[1])) + 360.0) % 360.0

    disagreement = abs((axis_bearing - wind_bearing_deg + 180.0) % 360.0 - 180.0)
    elongation = float(np.sqrt(eigvals.max() / max(eigvals.min(), 1e-9)))

    verdict = {
        "quantifiable": disagreement <= MAX_AXIS_WIND_DISAGREEMENT_DEG,
        "reason": None if disagreement <= MAX_AXIS_WIND_DISAGREEMENT_DEG else "axis_wind_mismatch",
        "pixels": pixels,
        "axis_bearing_deg": round(float(axis_bearing), 1),
        "wind_bearing_deg": round(wind_bearing_deg, 1),
        "axis_disagreement_deg": round(float(disagreement), 1),
        "elongation": round(elongation, 2),
        "length_m": round(float(np.sqrt(eigvals.max()) * 4.0 * pixel_size_m), 1),
    }
    log.info("ch4.preflight", **verdict)
    return verdict


def ime_threshold_plateau(enhancement: np.ndarray, mask_seed: np.ndarray,
                          pixel_area_m2: float, sigma: float) -> dict:
    """Sweep the mask threshold and look for a plateau.

    A plume that is genuinely separable from the background shows an IME that is
    flat over a range of thresholds. No plateau means the mask boundary is doing
    the arithmetic, not the plume.
    """
    results = []
    for k in (2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0):
        mask = ndimage.binary_dilation(mask_seed) & (enhancement > k * sigma)
        ime = float(enhancement[mask].sum()) * pixel_area_m2 * 16.04 * 1e-9 * 40.87
        results.append({"k_sigma": k, "ime_kg": round(ime, 2), "pixels": int(mask.sum())})

    values = np.array([r["ime_kg"] for r in results[1:-1]])   # ignore the extremes
    spread = float(values.ptp() / max(values.mean(), 1e-9))
    plateau = spread < 0.15

    log.info("ch4.ime.plateau", plateau=plateau, relative_spread=round(spread, 3),
             sweep=results)
    return {"plateau": plateau, "relative_spread": round(spread, 3),
            "ime_kg": round(float(values.mean()), 2), "sweep": results}

The plateau sweep is the single most informative diagnostic in the whole quantification chain, because it answers a question no single-threshold estimate can: is this number a property of the plume, or of the number I picked?

Deterministic Transformation Logic #

With the plume cleared for quantification, both estimators run over the same enhancement field. The IME implementation uses the plateau-averaged mass and an effective wind from an instrument-specific scaling; the CSF implementation integrates across transects perpendicular to the fitted plume axis and reports the median of the per-transect estimates along with their spread.

from dataclasses import dataclass

import numpy as np
import structlog

log = structlog.get_logger()

CH4_MOLAR_MASS = 16.04     # g mol-1
PPM_M_TO_KG_M2 = 1e-9 * 40.87 * CH4_MOLAR_MASS   # ppm-m -> kg m-2 at STP column


@dataclass(frozen=True)
class Flux:
    method: str
    q_kg_h: float
    q_sigma_kg_h: float
    u_eff_m_s: float
    detail: dict


def effective_wind(u10_m_s: float, instrument: str) -> tuple[float, float]:
    """Uэ = a*U10 + b, fitted per instrument against large-eddy simulations.

    The coefficients absorb each retrieval's spatial smoothing, so they are NOT
    transferable between instruments — using AVIRIS coefficients on Sentinel-2
    is a documented ~15% bias.
    """
    scalings = {
        "aviris-ng": (0.34, 0.44, 0.50),
        "emit": (0.33, 0.45, 0.55),
        "prisma": (0.30, 0.51, 0.60),
        "sentinel-2": (0.28, 0.60, 0.70),
    }
    if instrument not in scalings:
        raise ValueError(f"no effective-wind scaling for {instrument!r}; do not guess")
    a, b, rel_error = scalings[instrument]
    return a * u10_m_s + b, rel_error


def ime_flux(ime_kg: float, plume_length_m: float, u10_m_s: float,
             instrument: str, retrieval_rel_error: float = 0.20) -> Flux:
    """Q = IME * Ueff / L. The residence time L/Ueff is the transport model."""
    u_eff, wind_rel = effective_wind(u10_m_s, instrument)
    q = ime_kg * u_eff / plume_length_m * 3600.0

    # Wind dominates; combining in quadrature is the minimum honest treatment.
    rel = float(np.sqrt(retrieval_rel_error ** 2 + wind_rel ** 2))

    log.info("ch4.flux.ime", q_kg_h=round(q, 1), ime_kg=round(ime_kg, 2),
             length_m=round(plume_length_m, 1), u_eff=round(u_eff, 2),
             rel_uncertainty=round(rel, 3))
    return Flux("ime", round(q, 1), round(q * rel, 1), round(u_eff, 2),
                {"ime_kg": ime_kg, "plume_length_m": plume_length_m,
                 "retrieval_rel_error": retrieval_rel_error, "wind_rel_error": wind_rel})


def csf_flux(enhancement: np.ndarray, mask: np.ndarray, axis_bearing_deg: float,
             pixel_size_m: float, u10_m_s: float, instrument: str,
             transect_offsets_m: tuple[float, ...] = (120, 240, 360, 480)) -> Flux:
    """Integrate the column across transects perpendicular to the plume axis.

    Each transect is an independent estimate of the same source rate, so their
    spread is an empirical uncertainty that owes nothing to a stated error model.
    """
    u_eff, wind_rel = effective_wind(u10_m_s, instrument)

    ys, xs = np.nonzero(mask)
    weights = enhancement[ys, xs]
    cy, cx = np.average(ys, weights=weights), np.average(xs, weights=weights)

    theta = np.radians(axis_bearing_deg)
    along = np.array([np.sin(theta), -np.cos(theta)])     # unit vector down-plume
    across = np.array([along[1], -along[0]])

    estimates = []
    for offset in transect_offsets_m:
        step = offset / pixel_size_m
        centre = np.array([cx, cy]) + along * step
        # Sample the column along the perpendicular, half a pixel apart.
        samples = []
        for t in np.arange(-40.0, 40.0, 0.5):
            px, py = centre + across * t
            ix, iy = int(round(px)), int(round(py))
            if 0 <= iy < enhancement.shape[0] and 0 <= ix < enhancement.shape[1]:
                samples.append(enhancement[iy, ix])
        if len(samples) < 20:
            continue
        # Integrated column across the transect, in kg per metre of transect.
        line_density = float(np.sum(samples)) * PPM_M_TO_KG_M2 * (0.5 * pixel_size_m)
        estimates.append(line_density * u_eff * 3600.0)

    if not estimates:
        raise ValueError("no usable transects; plume too close to the scene edge")

    q = float(np.median(estimates))
    spread_rel = float(np.std(estimates, ddof=1) / max(q, 1e-9)) if len(estimates) > 1 else 0.3
    rel = float(np.sqrt(spread_rel ** 2 + wind_rel ** 2))

    log.info("ch4.flux.csf", q_kg_h=round(q, 1), transects=len(estimates),
             per_transect=[round(e, 1) for e in estimates],
             spread_rel=round(spread_rel, 3), rel_uncertainty=round(rel, 3))
    return Flux("csf", round(q, 1), round(q * rel, 1), round(u_eff, 2),
                {"per_transect_kg_h": [round(e, 1) for e in estimates],
                 "transect_spread_rel": round(spread_rel, 3), "wind_rel_error": wind_rel})


def reconcile(ime: Flux, csf: Flux, tolerance: float = 0.30) -> dict:
    """Report both, flag disagreement. Two structurally different estimators
    agreeing is stronger evidence than either one's stated interval."""
    mean = (ime.q_kg_h + csf.q_kg_h) / 2.0
    disagreement = abs(ime.q_kg_h - csf.q_kg_h) / max(mean, 1e-9)
    agreed = disagreement <= tolerance

    if not agreed:
        log.warning("ch4.flux.disagreement", ime_kg_h=ime.q_kg_h, csf_kg_h=csf.q_kg_h,
                    disagreement=round(disagreement, 3),
                    hint="ragged mask hurts IME; a bent axis hurts CSF")

    return {"q_kg_h": round(mean if agreed else min(ime.q_kg_h, csf.q_kg_h), 1),
            "q_sigma_kg_h": round(max(ime.q_sigma_kg_h, csf.q_sigma_kg_h), 1),
            "ime_kg_h": ime.q_kg_h, "csf_kg_h": csf.q_kg_h,
            "methods_agree": agreed, "disagreement": round(disagreement, 3),
            "u_eff_m_s": ime.u_eff_m_s, "reported": "mean" if agreed else "conservative_min"}

Note what reconcile does when the estimators disagree: it reports the lower of the two, not the mean, and marks the record. Conservativeness under uncertainty is a requirement of every carbon methodology, and it is far easier to defend a number chosen by a stated rule than one produced by averaging two figures you have just admitted are inconsistent.

Compliance Gating & Audit Trail Generation #

A flux estimate becomes reportable when it carries the evidence needed to re-derive it. Six fields are the minimum: the retrieval algorithm and version, the mask threshold and whether a plateau was found, the wind source with its timestamp and the effective-wind scaling used, both method estimates and their agreement flag, the uncertainty decomposition rather than a single combined number, and the steady-state verdict. Store them on the record, not only in the log, so a replay does not depend on log retention — the contract described in the MRV data schema reference.

Two gates then apply before the number reaches an inventory. First, conservativeness: where methods disagree beyond tolerance, or steady state is doubtful, report the lower estimate and mark the record, exactly as reconcile does. Second, materiality-scaled scrutiny: where a single plume contributes more than a few per cent of a facility’s reported total, the wind term should be replaced by an on-site measurement or the estimate should be repeated on a second acquisition before it is used. The uncertainty budget above shows why — the difference between a reanalysis wind and a mast is the difference between ±56% and ±27% on the same observation.

The annualisation step deserves its own record. Converting instantaneous fluxes to an annual mass requires a stated temporal model, and its assumptions belong in lineage alongside the measurements, wired through MRV data lineage and provenance tracking. An auditor who can see three observations and the persistence assumption can reconstruct the annual figure; one who sees only the annual figure cannot check anything.

Production Integration #

1. Ingest the enhancement field, plume mask, scene metadata, and wind record, refusing any scene whose CRS is missing — plume geometry drives the facility attribution and a wrong datum reassigns the emission.
2. Pre-flight with the checks above: reject unresolved plumes, plumes touching the scene edge, and masks whose principal axis disagrees with the wind bearing.
3. Sweep the threshold and take the plateau-averaged IME. If no plateau exists, record the plume as detected but unquantifiable rather than picking a threshold that produces a comfortable number.
4. Quantify twice — IME and CSF — using the instrument-specific effective-wind scaling, and reconcile with the conservative rule.
5. Decompose the uncertainty into wind, retrieval, mask, length, and steady-state terms, storing the decomposition rather than only its quadrature sum.
6. Emit one record per scene-facility pair, including non-detections with their detection limits, and hand it to the validation gates described under emissions data quality validation gates.

For batch operation, quantification is cheap relative to retrieval — the transects and threshold sweep are small array operations — so the sensible partition is per scene, with plumes processed in-process rather than fanned out. Cache the effective-wind scaling table as a versioned artefact, because changing a coefficient changes every historical flux and that is a restatement, not a bug fix.

Frequently Asked Questions #

Which method should I report if I can only report one? #

Report IME for compact plumes and CSF for long, well-developed ones — but the better answer is to compute both always and report the reconciled value with the agreement flag. The marginal cost is a few milliseconds; the marginal evidence is substantial. When forced to a single number by a submission template, use IME with the plateau check, since it degrades more gracefully at the small plume sizes that dominate most surveys.

Can I use the 10 m reanalysis wind directly instead of an effective-wind scaling? #

Only if you want a biased estimate. The 10 m wind is not the wind that transports a plume occupying the lowest tens of metres of a sheared boundary layer, and using it directly typically overestimates flux because Uₑ is usually below U₁₀ in the relevant range. The scaling coefficients are instrument-specific because they absorb the retrieval’s spatial smoothing; take them from a published fit for your instrument, record which fit you used, and never transfer coefficients between sensors.

What does it mean when IME and CSF disagree by a factor of two? #

Almost always a geometry problem rather than a physics one. A ragged or over-grown mask inflates IME because it sums noise over a larger area, while CSF is comparatively insensitive to the mask edge. A bent plume axis — common in shifting wind — breaks CSF because the transects are no longer perpendicular to the flow, while IME does not care about the axis. Check the elongation and axis-disagreement diagnostics first; if both look healthy, suspect the plume length used in IME.

How do I turn several instantaneous fluxes into an annual figure? #

With an explicit, stated temporal model. For a source known to be continuous, a mean of observations with an interval reflecting both measurement and temporal sampling error is defensible. For intermittent sources, model the duty cycle — from operational records where available, or from the detection frequency across observations, treating non-detections as censored observations below their detection limits. What is never defensible is multiplying one observed rate by 8760 hours without stating that you assumed persistence.

Does the plume mask threshold need to be the same across a survey? #

The sigma multiple should be consistent, but the absolute threshold will differ per scene because the noise floor differs. Fix the rule — for example, four times the scene’s robust noise estimate, with a plateau check — and apply it uniformly. Fixing the absolute ppm·m value instead makes your effective detection sensitivity vary with scene quality in a way that is invisible in the output and biases the survey towards clean scenes.

Related guides #

• Methane Plume Detection from Hyperspectral Imagery — the parent topic: retrieval, segmentation, and attribution.
• Troubleshooting False Methane Detections over Bright Surfaces — making sure the plume you are quantifying is real.
• Monte Carlo Uncertainty Propagation for Emission Factors — propagating this budget into a reported total.
• Emissions Data Quality & Validation Gates — the gates a flux record must pass to enter an inventory.