Biomass Estimation from LiDAR & SAR Fusion #

Biomass Estimation from LiDAR & SAR Fusion is the multi-sensor synthesis stage that converts vertical canopy structure and radar backscatter into spatially explicit aboveground biomass (AGB), the most computationally demanding step inside the Spatial Modeling & Carbon Stock Validation framework. It consumes cloud-masked, temporally aligned products from the upstream satellite imagery processing pipeline and depends on deterministic CRS alignment established in the foundational MRV layer, then hands calibrated AGB rasters and uncertainty bands downstream to ground-truth alignment and emission factor uncertainty mapping. Engineering a production-grade fusion stage means treating sensor disagreement, backscatter saturation, and sub-pixel drift as first-class, audit-traceable conditions rather than silent statistical noise.

Role in the MRV Workflow #

Within the Measurement, Reporting, and Verification (MRV) pipeline, fusion sits at the feature-synthesis boundary: it is downstream of ingestion and normalization, and upstream of calibration and compliance export. The two input modalities have complementary blind spots. Airborne or spaceborne LiDAR delivers precise vertical canopy structure but suffers from spatial gaps, high acquisition cost, and (for photon-counting spaceborne sensors) optical cloud dependency. Synthetic Aperture Radar (SAR) provides all-weather, wide-area coverage with sensitivity to woody volume and moisture, but experiences backscatter saturation in dense tropical canopies. Neither sensor alone produces a defensible wall-to-wall AGB surface, which is why fusion — not single-sensor extrapolation — is the only route to compliance-grade coverage.

The stage operates on co-registered raster tiles. LiDAR-derived Canopy Height Models (CHMs) typically align to a local datum or orthometric projection, whereas SAR backscatter composites are delivered in slant-range geometry or geocoded to a global reference grid (for example EPSG:4326 or EPSG:3857). Resolving that geometry into a single projected analysis CRS (commonly a local UTM zone such as EPSG:32633) is the precondition for every later operation, because any residual misalignment propagates directly into carbon accounting outputs.

Upstream dependencies are explicit: fusion expects terrain-corrected, multi-looked SAR and a hydrologically conditioned LiDAR ground model. Downstream consumers are equally explicit: the AGB raster and its per-pixel confidence interval feed plot-to-pixel calibration, after which threshold tuning for carbon stock baselines converts the validated surface into stable-versus-degraded masks for credit issuance. Because the stage straddles structural and radar sensing, its contract with the rest of the pipeline is unforgiving: it must emit not just a best estimate but a quantified, spatially resolved uncertainty layer that auditors can interrogate tile by tile.

Core Failure Modes #

Three failure modes account for the overwhelming majority of silent biomass errors in production fusion pipelines. Each has a concrete root cause and a measurable impact on reported carbon.

1. Sub-pixel co-registration drift. LiDAR CHMs and SAR backscatter rarely share a pixel grid exactly, and even a 0.3–0.5 pixel shift smears canopy edges against radar texture, decorrelating the height-to-backscatter relationship the model relies on. At 10 m SAR resolution this manifests as false biomass gradients along forest boundaries and systematic bias of 8–15% AGB in fragmented landscapes. Phase-correlation cross-matching followed by affine refinement, gated on a mutual-information score, is the only reliable defense; tiles that exceed the drift tolerance must be rejected, not “best-effort” stacked.

2. SAR backscatter saturation. C-band and L-band σ⁰ saturate logarithmically around 150–300 t/ha AGB, beyond which additional woody volume produces negligible change in radar reflectivity. A model that trusts SAR uniformly across the biomass range will compress high-density forest toward the saturation ceiling, understating stored carbon in exactly the tropical regions where credits are most valuable — observed underestimation routinely exceeds 30% in closed-canopy primary forest. The fix is regime-dependent weighting: SAR and texture features dominate the low-to-moderate regime, while LiDAR height percentiles take precedence above the saturation knee.

3. Uncoupled uncertainty propagation. Treating uncertainty as a global multiplier applied after prediction hides the fact that error structure is spatially heterogeneous: saturation zones, alignment-degraded tiles, and LiDAR-gap fallback regions each carry different variance. A flat ±15% band masks tiles that are genuinely outside tolerance and inflates confidence where it is not warranted, which is precisely the kind of defect third-party verifiers reject. Confidence intervals must be computed per pixel from bootstrapped ensemble variance and recomputed per jurisdictional boundary, never borrowed across regions.

Deterministic Implementation Architecture #

A scalable fusion stage requires distributed array processing, explicit CRS enforcement, deterministic fallback routing, and structured telemetry on every gate. The following Prefect flow leverages dask, xarray, rioxarray, and rasterio with structlog logging. Co-registration quality, saturation regime, and uncertainty are all surfaced as explicit, logged decisions rather than buried inside the model.

import structlog
import numpy as np
import rasterio
import rioxarray
import xarray as xr
import dask.array as da
from prefect import flow, task
from pyproj import CRS
from skimage.registration import phase_cross_correlation
from sklearn.ensemble import GradientBoostingRegressor

logger = structlog.get_logger()

TARGET_CRS = "EPSG:32633"        # local UTM analysis grid — area-true metres
DRIFT_TOL_PX = 0.5               # max admissible sub-pixel co-registration shift
SAT_KNEE_DB = -5.0               # σ⁰ (dB) proxy for ~250 t/ha SAR saturation
MAX_UNCERTAINTY_PCT = 15.0       # Verra VM0047 / 90% CI compliance ceiling


@task
def load_and_enforce_crs(lidar_path: str, sar_path: str) -> tuple[xr.DataArray, xr.DataArray]:
    """Open both modalities, reproject to a common projected CRS, fail loudly on mismatch."""
    chm = rioxarray.open_rasterio(lidar_path, chunks={"y": 2048, "x": 2048}).squeeze("band", drop=True)
    sigma0 = rioxarray.open_rasterio(sar_path, chunks={"y": 2048, "x": 2048}).squeeze("band", drop=True)

    if CRS.from_user_input(chm.rio.crs) != CRS.from_user_input(TARGET_CRS):
        logger.info("reproject_lidar", src=str(chm.rio.crs), dst=TARGET_CRS)
        chm = chm.rio.reproject(TARGET_CRS)
    if CRS.from_user_input(sigma0.rio.crs) != CRS.from_user_input(TARGET_CRS):
        logger.info("reproject_sar", src=str(sigma0.rio.crs), dst=TARGET_CRS)
        sigma0 = sigma0.rio.reproject(TARGET_CRS)

    # Snap SAR onto the LiDAR grid so every pixel is index-aligned downstream.
    sigma0 = sigma0.rio.reproject_match(chm)
    return chm, sigma0


@task
def coregistration_gate(chm: xr.DataArray, sigma0: xr.DataArray) -> float:
    """Phase-correlation drift check. Raises so Prefect can route to SAR fallback."""
    # Sample an overlapping, finite window to estimate residual sub-pixel shift.
    ref = np.nan_to_num(chm.isel(y=slice(0, 512), x=slice(0, 512)).values.astype("float32"))
    mov = np.nan_to_num(sigma0.isel(y=slice(0, 512), x=slice(0, 512)).values.astype("float32"))
    shift, error, _ = phase_cross_correlation(ref, mov, upsample_factor=20)
    drift_px = float(np.hypot(*shift))
    logger.info("coregistration_gate", drift_px=round(drift_px, 3), phase_error=round(float(error), 4))
    if drift_px > DRIFT_TOL_PX:
        logger.error("alignment_failed", drift_px=drift_px, tolerance=DRIFT_TOL_PX)
        raise ValueError(f"Sub-pixel drift {drift_px:.3f}px exceeds {DRIFT_TOL_PX}px; route to SAR fallback")
    return drift_px


@task
def build_fusion_features(chm: xr.DataArray, sigma0: xr.DataArray) -> xr.Dataset:
    """Regime-aware predictor matrix: weight LiDAR up where SAR saturates."""
    sigma0_db = 10.0 * xr.apply_ufunc(da.log10, sigma0.where(sigma0 > 0), dask="allowed")
    saturated = sigma0_db > SAT_KNEE_DB                 # high-biomass regime
    lidar_weight = xr.where(saturated, 0.8, 0.3)        # LiDAR leads past the knee
    sar_weight = 1.0 - lidar_weight

    features = xr.Dataset({
        "chm_p90": chm.rolling(y=5, x=5, center=True).reduce(np.nanpercentile, q=90),
        "chm_p75": chm.rolling(y=5, x=5, center=True).reduce(np.nanpercentile, q=75),
        "sigma0_db": sigma0_db,
        "lidar_weight": lidar_weight,
        "sar_weight": sar_weight,
    })
    logger.info("fusion_features", saturated_fraction=round(float(saturated.mean().compute()), 3))
    return features


@flow(name="biomass_estimation_lidar_sar_fusion", retries=2, retry_delay_seconds=30)
def run_fusion_pipeline(lidar_tile: str, sar_tile: str, model: GradientBoostingRegressor) -> xr.Dataset:
    logger.info("pipeline_start", lidar=lidar_tile, sar=sar_tile, target_crs=TARGET_CRS)

    chm, sigma0 = load_and_enforce_crs(lidar_tile, sar_tile)
    coregistration_gate(chm, sigma0)
    features = build_fusion_features(chm, sigma0)

    # Predict AGB per pixel; bootstrap the GBDT ensemble for a 90% interval.
    stack = xr.concat([features[v] for v in ("chm_p90", "chm_p75", "sigma0_db")], dim="feature")
    flat = stack.stack(px=("y", "x")).transpose("px", "feature").values
    preds = np.stack([est.predict(flat) for est in model.estimators_.ravel()])
    agb = preds.mean(axis=0).reshape(chm.shape)
    ci_lower = np.nanpercentile(preds, 5, axis=0).reshape(chm.shape)
    ci_upper = np.nanpercentile(preds, 95, axis=0).reshape(chm.shape)

    out = xr.Dataset(
        {
            "agb_t_ha": (("y", "x"), agb),
            "ci_lower_90": (("y", "x"), ci_lower),
            "ci_upper_90": (("y", "x"), ci_upper),
        },
        coords={"y": chm.y, "x": chm.x},
    )
    out["uncertainty_pct"] = (out["ci_upper_90"] - out["ci_lower_90"]) / out["agb_t_ha"] * 100.0
    out["audit_flag"] = out["uncertainty_pct"] > MAX_UNCERTAINTY_PCT
    out.rio.write_crs(TARGET_CRS, inplace=True)

    logger.info(
        "pipeline_complete",
        mean_agb=round(float(np.nanmean(agb)), 1),
        flagged_pct=round(float(out["audit_flag"].mean().compute()) * 100, 2),
    )
    return out

The flow is deterministic in three ways that matter for verification. CRS is declared once and enforced on both modalities with reproject_match, so pixels are index-aligned before any arithmetic. The co-registration gate raises rather than silently degrading, letting the orchestrator route alignment-failed tiles to a SAR-dominant fallback estimator. And the 90% interval is derived from the bootstrapped ensemble spread, not a constant multiplier — every pixel carries its own width.

Validation, Debugging & Compliance Mapping #

The pipeline emits three audit-ready artifacts per tile: a spatially explicit AGB raster, 90% confidence-interval bounds, and an uncertainty_pct mask whose audit_flag isolates tiles above the compliance ceiling. These map directly onto regulatory requirements. The IPCC 2006 Guidelines for National Greenhouse Gas Inventories Tier 3 tier mandates spatially resolved biomass models with documented error propagation; the per-pixel interval satisfies that documentation. For voluntary markets, Verra VM0047 and ART-TREES require project-level uncertainty to remain below 15% at 90% confidence, which is exactly what the audit_flag layer enforces by routing exceedance tiles to manual review or supplemental field sampling. The deterministic CRS handling and rejection gates provide the reproducibility evidence that ISO 14064-3 verifiers expect, and the flagged-fraction telemetry feeds the spatially explicit disclosure CSRD ESRS E1 anticipates.

Production fusion pipelines fail silently when geometric or statistical assumptions break, so each assumption needs an explicit safeguard:

Per-pixel feature-importance scores exported alongside the AGB raster let auditors confirm that SAR saturation zones are not artificially inflating carbon credits and that LiDAR-dominant regions cross-check against independent canopy-closure metrics. This is the interface that the detailed point-cloud-to-SAR fusion procedure builds on when it documents the full feature-stack provenance for registry submission.

Conclusion #

By enforcing explicit CRS declaration, a rejecting co-registration gate, regime-aware sensor weighting, and per-pixel bootstrapped uncertainty, Biomass Estimation from LiDAR & SAR Fusion becomes a repeatable, compliance-grade component of carbon accounting infrastructure rather than an opaque model call. The stage bridges high-resolution structural sensing and wide-area radar coverage, delivering spatially explicit estimates that withstand third-party verification. For the operational detail behind each gate — temporal harmonization windows, fallback routing logic, and registry-ready provenance — continue to the in-depth guide on fusing LiDAR point clouds with SAR for biomass estimation.

Related #

• Spatial Modeling & Carbon Stock Validation — parent framework and end-to-end pipeline overview
• Ground Truth Alignment for Carbon Models — plot-to-pixel calibration of the AGB surface this stage produces
• Emission Factor Uncertainty Mapping — propagating the confidence bands downstream
• Threshold Tuning for Carbon Stock Baselines — converting validated AGB into stable-versus-degraded masks
• Fusing LiDAR Point Clouds with SAR for Biomass Estimation — the detailed implementation procedure