Fusing LiDAR Point Clouds with SAR for Biomass Estimation #

Fusing discrete airborne LiDAR returns with continuous radar backscatter is the highest-leverage — and most error-prone — operation inside the Biomass Estimation from LiDAR & SAR Fusion stage of the Spatial Modeling & Carbon Stock Validation framework. Done correctly, it produces a wall-to-wall aboveground biomass (AGB) surface with defensible per-pixel uncertainty; done carelessly, it bakes silent spatial drift and backscatter saturation directly into issued carbon credits. This guide walks through the deterministic, audit-traceable procedure for merging a LiDAR canopy height model (CHM) with Sentinel-1 / ALOS-2 σ⁰ into a single regression-ready feature stack.

Because every downstream consumer treats the fused raster as ground truth, the operation depends on deterministic CRS alignment established at ingestion and on cloud-screened optical context from the upstream satellite imagery processing pipeline. It emits calibrated AGB and confidence bands that feed ground-truth alignment and emission factor uncertainty mapping further along the chain. The engineering intent here is reproducibility under MRV scrutiny — Verra VM0048, ISO 14064-3, and GHG Protocol spatial validation thresholds — not merely a statistically pleasing fit.

Root Cause Analysis #

The dominant failure mode in fusion pipelines is geometric, not statistical. LiDAR point clouds arrive in a local projected CRS (typically a UTM zone) with centimetre-level vertical accuracy, while SAR products are geocoded to EPSG:4326 with terrain-corrected radiometry. Any residual misalignment greater than roughly 0.5 pixels at 10 m SAR resolution decorrelates the height-to-backscatter relationship the model relies on, smearing canopy edges against radar texture and producing systematic AGB bias of 8–15% in fragmented landscapes. This is exactly the class of silent datum drift that auditors reject, because the error is invisible in summary statistics yet structurally present along every forest boundary.

The second root cause is temporal. SAR backscatter fluctuates with soil moisture and phenology, whereas LiDAR captures a single structural snapshot. Stacking a wet-season SAR scene against a dry-season LiDAR overflight injects moisture variance into a model that is supposed to read woody volume. Acquisitions must therefore be gated to a narrow window — commonly ±30 days — with scenes rejected when precipitation or soil-moisture anomalies exceed tolerance.

The third is physical: C-band and L-band σ⁰ saturate logarithmically around 150–300 t/ha, beyond which additional biomass produces negligible radar response. A model that trusts SAR uniformly across the biomass range compresses high-density forest toward the saturation ceiling, understating stored carbon by 30%+ in closed-canopy primary forest. The fix is regime-dependent weighting rather than a single global fit. The table below summarises the three causes and their detectable signatures.

Diagnostic Pipeline & Pre-Flight Validation #

Before any reprojection or resampling, the inputs must be inspected and any failure condition surfaced as an explicit, logged decision. The pre-flight gate below checks CRS validity, spatial overlap, LiDAR point density, and temporal proximity, emitting structlog telemetry on every check so the audit trail records why a tile was admitted or rejected rather than only the outcome.

import structlog
import rasterio
from pyproj import CRS, Transformer
from shapely.geometry import box
from datetime import datetime, timedelta

logger = structlog.get_logger()

TARGET_CRS = "EPSG:32618"        # LiDAR UTM analysis grid — area-true metres
MAX_DRIFT_PX = 0.5               # admissible sub-pixel co-registration shift
MIN_POINT_DENSITY = 5.0          # pts/m^2 floor for usable CHM
MAX_TEMPORAL_DELTA = timedelta(days=30)


def preflight_validate(lidar_meta: dict, sar_meta: dict) -> tuple[bool, dict]:
    """Detect fusion failure conditions before transformation. Returns (ok, report)."""
    report: dict = {"checks": {}, "rejected_reasons": []}

    # 1. CRS must be declared and parseable on both inputs — never assume.
    try:
        lidar_crs = CRS.from_user_input(lidar_meta["crs"])
        sar_crs = CRS.from_user_input(sar_meta["crs"])
    except Exception as exc:                       # noqa: BLE001
        logger.error("crs_unparseable", error=str(exc))
        return False, {"rejected_reasons": ["unparseable_crs"]}

    # 2. Spatial overlap in the shared analysis CRS (always_xy for lon/lat order).
    tx = Transformer.from_crs(sar_crs, lidar_crs, always_xy=True)
    sx0, sy0 = tx.transform(*sar_meta["bounds"][:2])
    sx1, sy1 = tx.transform(*sar_meta["bounds"][2:])
    overlap = box(*lidar_meta["bounds"]).intersection(box(sx0, sy0, sx1, sy1))
    frac = overlap.area / box(*lidar_meta["bounds"]).area if overlap.area else 0.0
    report["checks"]["overlap_fraction"] = round(frac, 4)
    if frac < 0.90:
        report["rejected_reasons"].append("insufficient_spatial_overlap")

    # 3. LiDAR point-density floor — sparse returns yield unreliable CHM percentiles.
    density = lidar_meta["point_count"] / max(overlap.area, 1.0)
    report["checks"]["point_density_pts_m2"] = round(density, 2)
    if density < MIN_POINT_DENSITY:
        report["rejected_reasons"].append("lidar_density_below_floor")

    # 4. Temporal gating — phenology/moisture divergence corrupts the σ⁰ signal.
    delta = abs(datetime.fromisoformat(sar_meta["datetime"])
                - datetime.fromisoformat(lidar_meta["datetime"]))
    report["checks"]["temporal_delta_days"] = delta.days
    if delta > MAX_TEMPORAL_DELTA:
        report["rejected_reasons"].append("temporal_window_exceeded")

    ok = not report["rejected_reasons"]
    logger.info("preflight_complete", ok=ok, **report["checks"],
                reasons=report["rejected_reasons"])
    return ok, report

Temporal gating should be enforced at metadata ingestion using STAC datetime query bounds and an external meteorological API, so non-compliant scenes are dropped before rasterization rather than after expensive resampling.

Deterministic Transformation Logic #

Once a tile passes pre-flight, the SAR layer is orthorectified onto the LiDAR analysis grid and the two modalities are reduced to orthogonal structural and dielectric features. Determinism is non-negotiable: the same inputs must always yield byte-identical outputs, which means pinned resampling kernels, explicit always_xy axis ordering, and area-preserving percentile aggregation rather than naive mean downsampling.

import numpy as np
import rasterio
import rioxarray
import xarray as xr
from rasterio.warp import calculate_default_transform, reproject, Resampling
from scipy.ndimage import uniform_filter

logger = structlog.get_logger()


def fuse_lidar_sar(sar_path: str, chm_1m_path: str, out_path: str,
                   target_crs: str = TARGET_CRS, res_m: float = 10.0) -> str:
    """Orthorectify SAR to the LiDAR grid and assemble a regression-ready stack."""
    # 1. Reproject SAR onto the LiDAR projected CRS — cubic kernel, explicit transform.
    with rasterio.open(sar_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, target_crs, src.width, src.height, *src.bounds, resolution=res_m)
        dst_meta = {**src.meta, "crs": target_crs, "transform": transform,
                    "width": width, "height": height}
        sar = np.empty((src.count, height, width), dtype="float32")
        for i in range(1, src.count + 1):
            reproject(source=rasterio.band(src, i), destination=sar[i - 1],
                      src_transform=src.transform, src_crs=src.crs,
                      dst_transform=transform, dst_crs=target_crs,
                      resampling=Resampling.cubic)
    logger.info("sar_reprojected", target_crs=target_crs, shape=sar.shape)

    # 2. CHM 1 m -> 10 m using the 75th percentile to preserve upper-canopy structure.
    chm_10m = (rioxarray.open_rasterio(chm_1m_path)
               .rio.reproject(target_crs, resolution=res_m,
                              resampling=Resampling.q3))

    # 3. Refined-Lee speckle suppression on the VH intensity band (index 1).
    vh = sar[1]
    mean = uniform_filter(vh, size=5)
    var = uniform_filter(vh ** 2, size=5) - mean ** 2
    noise_var = 0.09                                # Sentinel-1 GRD nominal ENL
    vh_filtered = mean + (var / (var + noise_var)) * (vh - mean)

    # 4. Log-transform the cross-pol ratio to stabilise variance in dense canopy.
    vv = np.clip(sar[0], 1e-6, None)
    vh_vv = np.log(np.clip(vh_filtered, 1e-6, None) / vv)

    # 5. Assemble the feature stack on a single, shared grid.
    coords = {"y": chm_10m.y, "x": chm_10m.x}
    stack = xr.Dataset({
        "chm":   (["y", "x"], chm_10m.squeeze().values),
        "vh":    (["y", "x"], vh_filtered),
        "vh_vv": (["y", "x"], vh_vv),
    }, coords=coords)
    stack.rio.write_crs(target_crs, inplace=True)
    stack.to_netcdf(out_path)
    logger.info("fusion_stack_written", path=out_path, vars=list(stack.data_vars))
    return out_path

The resulting multi-dimensional array is the input to a regime-weighted allometric regression or gradient-boosted ensemble. Keeping the cross-polarization ratio log-transformed before stacking is what stabilises variance across the saturation knee.

Compliance Gating & Audit Trail Generation #

MRV compliance demands explicit, spatially resolved uncertainty. Biomass estimates must carry per-pixel bounds satisfying ISO 14064-3 verification thresholds — typically ≤10% uncertainty at 90% confidence for project-scale baselines. Uncertainty is propagated through the allometric form AGB = exp(β₀ + β₁·ln(CHM) + β₂·ln(VH/VV) + ε) using the law of propagation of uncertainty, u_c² = Σ(∂AGB/∂x_i)²·u(x_i)² + 2·ΣΣ(∂AGB/∂x_i)(∂AGB/∂x_j)·cov(x_i,x_j), evaluated by first-order Taylor expansion or Monte Carlo sampling.

The gate halts execution when residual spatial autocorrelation (Moran’s I) exceeds 0.35, when VH σ⁰ breaches the tropical saturation threshold of -12 dB, or when LiDAR density falls below 5 pts/m² across more than 15% of the area. Every run serialises its inputs, transforms, and outcomes into an append-only audit record with SHA-256 checksums, mapped to the MRV data lineage requirements that the carbon credit registry submission depends on.

import json
import hashlib
from datetime import datetime

logger = structlog.get_logger()


def _sha256(path: str) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(1 << 20), b""):
            h.update(chunk)
    return h.hexdigest()


def gate_and_audit(run_id: str, uncertainty_raster, inputs: dict,
                   outputs: dict, params: dict, threshold_pct: float = 10.0) -> dict:
    """Enforce the compliance ceiling and emit a signed, append-only audit record."""
    mean_u = float(uncertainty_raster.mean())
    passed = mean_u <= threshold_pct

    audit = {
        "pipeline_version": "v2.4.1-mrv",
        "run_id": run_id,
        "execution_timestamp": datetime.utcnow().isoformat() + "Z",
        "spatial_reference": {"target_crs": TARGET_CRS,
                              "dem_version": "Copernicus_30m_v2023"},
        "inputs": {k: {"stac_id": v["id"], "sha256": _sha256(v["path"])}
                   for k, v in inputs.items()},
        "parameters": params,
        "compliance_gating": {
            "mean_pixel_uncertainty_pct": round(mean_u, 3),
            "passes_threshold": passed,
            "iso_14064_3_compliant": mean_u <= 10.0,
            "verra_vm0048_ceiling_pct": threshold_pct,
        },
        "outputs": {k: {"path": v, "sha256": _sha256(v)} for k, v in outputs.items()},
    }
    with open(f"audit_{run_id}.json", "w") as fh:
        json.dump(audit, fh, indent=2)

    if not passed:
        logger.error("compliance_gate_failed", run_id=run_id, mean_u=mean_u)
        raise RuntimeError(f"Uncertainty {mean_u:.2f}% exceeds {threshold_pct}% ceiling")
    logger.info("compliance_gate_passed", run_id=run_id, mean_u=mean_u)
    return audit

Any deviation from default allometric coefficients requires version-controlled justification and auditor sign-off before deployment, and gap-interpolation regions (kriging or random-forest imputation) must be flagged with explicit coverage percentages so verifiers can isolate modelled-versus-observed pixels.

Production Integration #

A production fusion run is a strictly ordered, containerised sequence — Docker or Singularity with pinned dependency versions — so that every estimated tonne is reproducible and audit-ready. Batch execution chunks tiles through a dask-backed scheduler, but each tile follows the same six-step contract:

1. Ingest — pull SAR and LiDAR products by STAC ID, verify SHA-256 checksums, and record source provenance.
2. Diagnose — run preflight_validate; reject and log any tile that fails CRS, overlap, density, or temporal gates before spending compute on resampling.
3. Transform — execute fuse_lidar_sar to orthorectify SAR onto the LiDAR grid and assemble the feature stack with deterministic kernels.
4. Validate — predict AGB with regime weighting, propagate per-pixel uncertainty, and evaluate Moran’s I and saturation diagnostics.
5. Export — write the AGB raster and confidence band as cloud-optimised GeoTIFFs to versioned, immutable object storage.
6. Submit — call gate_and_audit; on pass, attach the signed audit JSON and forward to the registry submission queue, threshold tuning handled downstream by threshold tuning for carbon stock baselines.

Chunked I/O matters at scale: read SAR scenes windowed to the LiDAR tile footprint rather than whole-scene, and stream NetCDF stacks lazily through xarray so memory stays bounded across thousands of tiles. With this contract in place, every fused biomass surface is traceable from raw return to issued credit.

Related #

• Biomass Estimation from LiDAR & SAR Fusion — parent topic and stage overview
• Ground-Truth Alignment for Carbon Models — plot-to-pixel calibration of fused outputs
• Emission Factor Uncertainty Mapping — propagating fusion uncertainty into emission factors
• Threshold Tuning for Carbon Stock Baselines — converting validated AGB into stable/degraded masks
• Geospatial CRS Alignment — the ingestion-stage prerequisite for fusion