Spatial Pipeline EngineeringMRV Automation

Monthly Temporal Aggregation of NDVI for Land Cover Change

In automated MRV (Measurement, Reporting, and Verification) architectures, raw daily reflectance observations are statistically unusable for regulatory carbon accounting. This guide is the task-level recipe under Temporal Aggregation for Land-Use Change, the compositing discipline within the Satellite Imagery Processing for Emissions Tracking stack. It depends on the clear-sky masks produced upstream by Sentinel-2 and Landsat cloud masking workflows and feeds the rolling baselines consumed downstream by deforestation alert generation pipelines.

Atmospheric interference, sensor geometry drift, and acquisition gaps introduce stochastic noise that obscures structural vegetation transitions. Monthly temporal aggregation of NDVI for land cover change resolves this by transforming heterogeneous pixel-level observations into deterministic, compliance-ready composites. The engineering objective is not statistical smoothing; it is the construction of an auditable temporal baseline that survives ESG scrutiny, aligns with GHG Protocol activity data requirements, and scales across multi-sensor tile streams as activity data entering the wider MRV architecture and carbon accounting fundamentals stack.

Monthly NDVI aggregation flow with observation-density branch A daily NDVI stack from Sentinel-2 and Landsat passes through a QA clear-sky filter keyed on SCL or QA_PIXEL, then a pinned monthly resample using a median or 75th-percentile reducer. A decision gate asks whether each cell has at least three clear observations per month: cells that pass become a deterministic monthly composite; cells that fail route to a 90-day rolling median fallback with the interpolation explicitly flagged. Both paths converge into a single audit manifest recording QA pass rate, datum and lineage. Daily NDVI stack QA clear-sky filter Monthly resample Monthly composite Rolling median fallback Audit manifest Sentinel-2 · Landsat SCL / QA_PIXEL median / p75 · 1MS deterministic baseline 90-day · interpolation flagged QA rate · datum · lineage ISO 14064-3 ready ≥ 3 clear obs per month? yes no

Root Cause Analysis: Why Daily Stacks Fail Verification

A naive pipeline that differences raw daily NDVI scenes, or that means them without QA gating, cannot satisfy third-party verification. Three structural problems drive the failure.

First, residual atmospheric contamination is directional, not random. Cloud, shadow, and thin cirrus depress near-infrared reflectance in exactly the way real canopy loss does, so the error does not average out — it biases the composite toward false disturbance. Even 15% residual cirrus contamination can depress NDVI by 0.08–0.12, enough to trip a deforestation threshold on a pixel that never changed. This is why a decoded clear-sky mask, the contract established upstream by the cloud masking workflows, is a precondition rather than a refinement.

Second, irregular revisit creates unequal evidentiary weight. Tropical tiles may yield two clear acquisitions in a wet month and twenty in a dry one. A reducer that ignores observation density treats a single haze-fringed pixel with the same authority as a twenty-deep stack, and a verifier cannot distinguish a confident measurement from a lucky guess. Without an explicit minimum-observation gate the composite silently encodes that imbalance.

Third, non-deterministic reduction breaks reproducibility. Floating window boundaries, unpinned reducers, or NDVI_max compositing (which selects the single brightest — often most contaminated — observation) produce outputs that cannot be regenerated identically. ISO 14064-3 verification demands that the same inputs and parameters yield the same composite, every run. The fix is a pinned calendar-month window, a robust median or 75th-percentile reducer that suppresses transient soil-moisture and BRDF edge effects, and an embedded record of every parameter used.

Diagnostic Pipeline: Pre-Flight Observation Sufficiency

Before any reduction runs, inspect the inputs and detect the failure conditions above. The pre-flight gate confirms a machine-readable CRS, a monotonic time axis, and — critically — that each calendar month carries enough clear-sky observations to support a defensible composite. Strict CRS alignment is validated here because a datum mismatch relocates the entire reporting cell and a sub-pixel offset between epochs is indistinguishable from real surface change once observations are reduced into one composite.

import xarray as xr
import numpy as np
import pandas as pd
import structlog

log = structlog.get_logger()

# Sentinel-2 SCL clear-sky classes: Vegetation, Bare Soil, Water, Unclassified
S2_SCL_CLEAR_CLASSES = {4, 5, 6, 7}
# Landsat C2 QA_PIXEL: clear requires cloud bit (6) and cloud-confidence bits (8-9) unset
LANDSAT_CLEAR_MASK = 0x0040 | 0x0300


def preflight_ndvi_stack(
    ndvi_stack: xr.DataArray,
    qa_mask: xr.DataArray,
    sensor: str = "S2",
    min_obs: int = 3,
    target_crs: str = "EPSG:4326",
) -> dict:
    """Inspect a daily NDVI stack and its QA band before monthly aggregation.

    Detects: missing/mismatched CRS, non-monotonic timestamps, and calendar
    months whose clear-sky observation count falls below the compliance floor.
    """
    issues: list[str] = []

    # 1. CRS must be declared and match the canonical target
    crs = getattr(ndvi_stack.rio, "crs", None)
    if crs is None:
        issues.append("missing_crs")
    elif str(crs) != target_crs:
        issues.append(f"crs_mismatch:{crs}!={target_crs}")

    # 2. Time axis must be strictly increasing (no duplicate acquisitions)
    times = pd.DatetimeIndex(ndvi_stack["time"].values)
    if not times.is_monotonic_increasing or times.has_duplicates:
        issues.append("non_monotonic_time_axis")

    # 3. Per-month clear-sky observation density
    if sensor == "S2":
        clear = xr.zeros_like(qa_mask, dtype=bool)
        for cls in S2_SCL_CLEAR_CLASSES:
            clear = clear | (qa_mask == cls)
    elif sensor in ("L8", "L9"):
        clear = (qa_mask & LANDSAT_CLEAR_MASK) == 0
    else:
        raise ValueError(f"Unsupported sensor: {sensor}. Use 'S2', 'L8', or 'L9'.")

    monthly_obs = clear.resample(time="1MS").sum(dim="time")
    thin_months = int((monthly_obs.min(dim=["y", "x"]) < min_obs).sum())

    report = {
        "sensor": sensor,
        "n_acquisitions": int(times.size),
        "months_below_min_obs": thin_months,
        "min_obs_threshold": min_obs,
        "sufficient": not issues and thin_months == 0,
        "issues": issues,
    }
    log.info("ndvi.preflight", **report)
    return report

A tile that reports sufficient=False is not discarded; the flag routes thin months to the rolling-fallback path in the transformation step and records why, so the gap is visible to a verifier rather than silently interpolated.

Observation count per monthly composite, and what it does to the interval A chart of the 95 percent confidence half-width on a monthly mean vegetation index against the number of clear observations contributing to it. With one observation the half-width is 0.11, with three it is 0.06, with six 0.04, with twelve 0.03. A shaded band marks the region below three observations where the interval exceeds the seasonal signal being measured. A panel notes that a composite carries its observation count as a band so that downstream aggregation can weight by it, and that treating a one-observation cell and a twelve-observation cell as equivalent hides a fourfold difference in precision. A composite is not one number, it is a number and a count 95% half-width on a monthly mean index, by contributing observations. interval > signal 0.12 0.08 0.04 1 3 6 12 clear observations in the month Carry n_obs as a band so downstream aggregation can weight by it. A one-observation cell and a twelve-observation cell differ fourfold in precision.

Deterministic Transformation Logic

The aggregation engine operates on chunked Cloud-Optimized GeoTIFFs using xarray and dask.array, enabling out-of-core processing across continental extents. While NDVI_max composites are standard for phenological tracking, land cover change detection for carbon accounting requires NDVI_median or the 75th percentile to suppress transient soil-moisture spikes and BRDF-induced edge effects. Multi-sensor harmonization is a prerequisite: Sentinel-2 MSI and Landsat 8/9 OLI must be normalized to a common reflectance scale using sensor-specific gain/offset tables before stacking, or an artificial NDVI step-change appears at every sensor handoff.

The routine below decodes the QA band, applies the clear-sky mask, reduces each pinned calendar month, and masks any cell that falls below the observation floor so a downstream 90-day rolling median can fill it explicitly.

import xarray as xr
import numpy as np
import dask.array as da
import pandas as pd
import rioxarray  # registers the xarray ".rio" accessor
import structlog
from datetime import datetime, timezone

log = structlog.get_logger()

# Landsat Collection 2 QA_PIXEL: clear requires bit 6 (cloud) unset
# and bits 8-9 (cloud confidence) both zero.
LANDSAT_CLEAR_MASK = 0x0040 | 0x0300  # bits to check for cloud contamination

# Sentinel-2 SCL valid (clear-sky) classes
S2_SCL_CLEAR_CLASSES = {4, 5, 6, 7}  # Vegetation, Bare Soil, Water, Unclassified


def aggregate_monthly_ndvi(
    ndvi_stack: xr.DataArray,
    qa_mask: xr.DataArray,
    sensor: str = "S2",
    method: str = "median",
    min_obs: int = 3,
    target_crs: str = "EPSG:4326",
) -> tuple[xr.DataArray, dict]:
    """
    Aggregate daily NDVI into monthly composites with QA-aware filtering,
    CRS validation, and compliance observation thresholds.

    qa_mask semantics:
      - Sentinel-2: integer SCL band (4,5,6,7 = clear)
      - Landsat C2:  integer QA_PIXEL band (clear when cloud bits are unset)

    Returns: (monthly_composite, audit_manifest)
    """
    # Distortion gate: enforce strict CRS alignment before any temporal stacking
    if ndvi_stack.rio.crs is not None and str(ndvi_stack.rio.crs) != target_crs:
        raise ValueError(f"CRS mismatch: {ndvi_stack.rio.crs} != {target_crs}")

    # Build boolean clear-sky mask from sensor-specific QA
    if sensor == "S2":
        # True where SCL class is in the valid set
        clear_mask = xr.zeros_like(qa_mask, dtype=bool)
        for cls in S2_SCL_CLEAR_CLASSES:
            clear_mask = clear_mask | (qa_mask == cls)
    elif sensor in ("L8", "L9"):
        # True where cloud bit (6) and cloud-confidence bits (8-9) are all zero
        clear_mask = (qa_mask & LANDSAT_CLEAR_MASK) == 0
    else:
        raise ValueError(f"Unsupported sensor: {sensor}. Use 'S2', 'L8', or 'L9'.")

    ndvi_clean = ndvi_stack.where(clear_mask)

    # Pinned calendar-month window ("1MS") + robust reducer for reproducibility
    if method == "max":
        monthly_composite = ndvi_clean.resample(time="1MS").max(dim="time")
    elif method == "p75":
        monthly_composite = ndvi_clean.resample(time="1MS").quantile(0.75, dim="time")
    else:
        monthly_composite = ndvi_clean.resample(time="1MS").median(dim="time")

    # Audit gate: enforce minimum clear-sky observation count per cell
    obs_count = clear_mask.resample(time="1MS").sum(dim="time")
    low_density_mask = obs_count < min_obs

    # Mask cells below threshold — downstream fallback applies rolling 90-day composite
    monthly_composite = monthly_composite.where(~low_density_mask)

    # Generate audit manifest for MRV verification
    audit_manifest = {
        "pipeline_version": "2.4.1-mrv",
        "sensor": sensor,
        "aggregation_method": method,
        "min_obs_threshold": min_obs,
        "target_crs": target_crs,
        "temporal_range": f"{ndvi_stack.time.values[0]} to {ndvi_stack.time.values[-1]}",
        "low_density_cells_pct": float((low_density_mask.sum() / low_density_mask.size) * 100),
        "qa_pass_rate_pct": float((clear_mask.sum() / clear_mask.size) * 100),
        "generated_utc": datetime.now(timezone.utc).isoformat(),
        "compliance_standard": "GHG Protocol Scope 3 LULUCF Activity Data",
    }

    # Stamp provenance attributes directly into the xarray dataset
    monthly_composite.attrs.update({
        "audit_manifest": audit_manifest,
        "qa_source": "SCL clear-sky classes {4,5,6,7} or Landsat QA_PIXEL cloud-bit check",
        "interpolation_flag": "90-day rolling median fallback applied where obs < min_obs",
    })

    log.info("ndvi.aggregated", **audit_manifest)
    return monthly_composite, audit_manifest

When acquisition density falls below the compliance threshold (typically <3 clear observations per month), fallback routing triggers a rolling 90-day median composite rather than propagating null values. This preserves temporal continuity for downstream change detection while explicitly flagging interpolated cells in the provenance metadata.

Compliance Gating & Audit Trail Generation

Regulatory carbon accounting requires transparent activity data lineage. The routine above embeds a machine-readable audit_manifest directly into the output xarray attributes. This manifest survives downstream serialization to Parquet or GeoTIFF and satisfies third-party verifier requirements under ISO 14064-3 and the GHG Protocol Corporate Accounting and Reporting Standard. Those attributes flow directly into MRV data lineage and provenance tracking, where they become the queryable record an auditor traces.

Key compliance gates enforced:

  1. Observation density threshold. Cells with <3 clear-sky acquisitions are masked and logged. Verifiers can trace interpolation flags to specific temporal windows rather than discovering silent gaps.
  2. QA pass rate tracking. The qa_pass_rate_pct metric documents atmospheric clearance performance per tile. Drops below 60% trigger automatic pipeline alerts for manual review.
  3. CRS and geolocation integrity. Strict projection enforcement prevents sub-pixel drift during temporal stacking, ensuring spatial consistency across multi-year baselines.
  4. Fallback transparency. When rolling 90-day composites are invoked, the interpolation_flag attribute explicitly marks affected pixels, preventing false attribution of carbon stock changes to unverified data.

For authoritative QA parsing specifications, consult the USGS Landsat Collection 2 Quality Assessment Band and the ESA Sentinel-2 Level-2A Scene Classification Map.

Production Integration

Deploy the routine within an async tile-processing framework, following a fixed ingest → diagnose → transform → validate → export → submit sequence:

  1. Ingest. Query the STAC API for Sentinel-2 MSI and Landsat 8/9 OLI items over the tile footprint and reporting window, reading chunked COGs lazily so out-of-core array operations never block on a single oversized stack. Throughput at continental scale leans on async satellite tile processing with Dask.
  2. Diagnose. Run preflight_ndvi_stack to confirm CRS, a monotonic time axis, and per-month observation sufficiency; route thin months to the fallback path.
  3. Transform. Call aggregate_monthly_ndvi with method="median" (or "p75" for soil-moisture-prone semi-arid tiles), reducing each pinned calendar month.
  4. Validate. Assert the qa_pass_rate_pct and low_density_cells_pct gates and reject any composite that breaches the configured floors.
  5. Export. Serialize to Parquet or GeoTIFF with the audit_manifest and interpolation_flag attributes intact.
  6. Submit. Feed the composites into a CUSUM or Bayesian change-detection module — the change models behind the deforestation alert generation pipelines — and forward the lineage to registry submission. The aggregated trajectories also feed parcel-level estimates in the spatial modeling and carbon stock validation layer.

Use dask.distributed to schedule chunked COG reads across distributed workers, routing STAC item lists through a priority queue that balances sensor availability and cloud-cover forecasts. By enforcing strict QA filtering, statistical robustness, and embedded audit trails, monthly temporal aggregation of NDVI transforms noisy optical observations into verifiable carbon accounting inputs that survive regulatory scrutiny and provide the deterministic foundation required for automated MRV compliance.

Where NDVI saturates, and which index to use instead A curve of index value against increasing leaf area index from 0 to 8. NDVI rises steeply to about 0.85 by a leaf area index of 3 and is then almost flat, so it cannot distinguish moderately dense canopy from very dense canopy — exactly the range where forest carbon varies most. EVI continues rising to a leaf area index of about 6 before flattening. A near-infrared reflectance of vegetation index continues further still. A shaded band marks the saturation region above a leaf area index of 3. A panel notes that NDVI remains excellent for detecting the presence or absence of vegetation and poor for grading its density, which is why change detection and biomass estimation should not share an index by default. NDVI saturates exactly where forest carbon varies Index response against leaf area index. NDVI saturation region high low 0 3 5.5 8 leaf area index NDVI EVI NIRv NDVI: presence or absence, yes. Grading density, no. Detection and biomass should not share an index.

Frequently Asked Questions

Why does the observation count matter as much as the composite value?

Because it is the composite’s precision, and precision varies enormously across a map. A monthly mean from one observation carries roughly four times the interval of one from twelve, and treating the two as equivalent in a downstream aggregation understates uncertainty exactly where the data was thinnest. Carry the count as a band, weight by it when aggregating, and report the distribution — a period whose median cell had two observations supports a much weaker claim than one whose median cell had ten.

Should NDVI be the default index for land-use change?

For detecting presence and absence of vegetation, yes; for grading density, no. NDVI saturates above a leaf area index of roughly three, which is precisely the range where forest carbon varies most, so it registers clearing sharply and degradation barely at all. Use it for change detection and use EVI, NIRv, or a structural signal for anything that needs to distinguish dense canopy from denser canopy.

How should monthly composites handle a month with no clear observations?

Mark the cell unobserved and let downstream decide. Interpolating across the gap is often necessary for a seasonal model and must be flagged as interpolated; substituting the previous month silently is not acceptable because it creates artificial persistence that suppresses real change. The fraction of cells requiring interpolation in a period is a headline quality statistic and belongs in the period’s metadata.

Does the compositing window need to be exactly a calendar month?

Not necessarily, but it must be consistent and documented. Calendar months are convenient for reporting alignment and awkward for phenology, since a growing season rarely respects month boundaries. Some pipelines use fixed 16-day or 32-day windows for analysis and aggregate to calendar months only for reporting, which keeps phenological comparability and reporting alignment separate rather than compromising both.

How do I compare composites across sensors?

Harmonise before compositing, never after. Sentinel-2 and Landsat differ in band centres, widths, and radiometric calibration, so their indices differ systematically over identical ground. Apply published cross-sensor transformation coefficients to bring one onto the other’s scale, validate the harmonisation over a stable reference area, and record which sensor contributed to each cell — a composite that silently mixes sensors carries a step change at every point where the mix shifts.

Should the composite store the contributing scene identifiers?

Yes, as an array column or a sidecar table keyed on the cell. Knowing which scenes produced a composite is what makes it reproducible and what lets a suspicious cell be traced back to a specific acquisition — usually the one whose mask leaked. The storage cost is small relative to the imagery, and the alternative is re-running the scene query later, which returns a different set once the archive has been reprocessed.

Should composites be recomputed when the archive is reprocessed?

Only deliberately, and as a restatement. A provider reprocessing its archive changes the inputs to every historical composite, so re-running produces different numbers for periods already reported. Keep the original composites and their pinned scene lists, decide explicitly whether to adopt the reprocessed archive, and if you do, recompute the whole affected history at once so the series stays internally consistent rather than changing at an arbitrary point.