Troubleshooting LiDAR–SAR Coregistration Drift #

Coregistration drift is the quietest way to poison a fused biomass product. When LiDAR-derived canopy height and SAR backscatter no longer sit on the same physical footprint, the height-to-backscatter relationship the model was trained on decorrelates, and aboveground biomass (AGB) error inflates without any single pixel looking obviously wrong. This guide is the corrective procedure under Biomass Estimation from LiDAR–SAR Fusion, the multi-sensor synthesis stage of the Spatial Modeling & Carbon Stock Validation framework. It exists because drift is not a modelling problem you can regularise away — it is a geometric fault that must be measured and reversed before any regression runs.

Left unaddressed, a half-pixel offset at 10 m resolution smears every canopy edge against radar texture, and a two-pixel offset can double the fused biomass RMSE across fragmented landscapes while summary statistics stay flat. The corrective loop below belongs upstream of ground-truth alignment for carbon models, which cannot calibrate plot-to-pixel agreement on a stack that is still spatially inconsistent. The engineering objective is a validated, coregistered feature stack whose residual shift is provably within a declared tolerance and whose correction is recorded as an auditable transform, not an ad-hoc nudge.

Root Cause Analysis #

Drift is a symptom with several distinct causes, and treating them as one leads to over-correction that manufactures its own error. The most common source is geolocation and orbit error in the SAR product. Sentinel-1 and ALOS-2 geocoding relies on a precise orbit ephemeris and a reference DEM; when the orbit state vectors are predicted rather than restituted, or when the scene is processed against a coarse global DEM, the geocoded pixel can land tens of metres from its true ground position. That offset is spatially coherent, so it presents as a clean translational shift rather than random scatter — which is precisely why cross-correlation can recover it.

Terrain effects are the second cause and the hardest to fully reverse. In steep or mountainous terrain, SAR layover and foreshortening displace backscatter towards the sensor, and the magnitude of that displacement varies with local slope and aspect. Unlike a global orbit offset, terrain-induced drift is spatially variable: a single (dy, dx) translation will reduce the mean error but leave residual, slope-correlated misalignment along ridgelines. When the terrain component dominates, the correct response is to re-run range-Doppler terrain correction against a better DEM, not to force a global shift that trades ridge error for valley error.

The third cluster of causes is grid and datum handling inside the pipeline itself. A DEM error propagates directly into the SAR geocoding solution; a datum mismatch — for example treating an EPSG:4326 product as if it were already on the LiDAR datum — introduces a systematic offset that looks like drift but is actually a silent transform bug, the same class of fault covered in debugging silent datum shifts in carbon pipelines. Naive resampling adds a fourth: nearest-neighbour reprojection onto a mismatched grid origin can shear the two layers by up to half a pixel before any true drift is even considered. Diagnosis therefore has to separate a coherent translational offset (correctable by a shift) from a spatially variable geometric distortion (not correctable by a shift) before deciding what to do.

Diagnostic Pipeline / Pre-Flight Validation #

The diagnostic step measures drift as a signed, subpixel quantity rather than asserting alignment on faith. It first enforces CRS and resolution parity — cross-correlation is only meaningful when both rasters occupy the same projection and pixel size — then reads the shared overlap window and estimates the shift with phase_cross_correlation, which locates the correlation peak in the frequency domain and refines it to subpixel precision. The returned (dy, dx) is the translation that aligns the moving SAR layer onto the reference CHM; its magnitude, compared against tolerance, decides whether the tile is admitted, corrected, or escalated.

import numpy as np
import rasterio
from rasterio.windows import from_bounds
from pyproj import CRS
from skimage.registration import phase_cross_correlation
import structlog

log = structlog.get_logger()

DRIFT_TOLERANCE_PX = 0.5      # max admissible residual coregistration shift
UPSAMPLE_FACTOR = 20          # subpixel resolution of the cross-correlation peak


def estimate_coregistration_drift(
    chm_path: str, sar_path: str, expected_res_m: float = 10.0
) -> dict:
    """Estimate the (dy, dx) pixel shift between a LiDAR CHM and SAR backscatter.

    Both rasters must already share a CRS and grid resolution. Returns a report
    carrying the subpixel offset in pixels and metres and a boolean flag when the
    shift magnitude exceeds the declared tolerance.
    """
    with rasterio.open(chm_path) as chm, rasterio.open(sar_path) as sar:
        # 1. CRS parity — a silent datum mismatch masquerades as physical drift.
        chm_crs, sar_crs = CRS.from_user_input(chm.crs), CRS.from_user_input(sar.crs)
        if chm_crs != sar_crs:
            log.error("crs_mismatch", chm=chm_crs.to_epsg(), sar=sar_crs.to_epsg())
            raise ValueError("Inputs must share a CRS before drift estimation")

        # 2. Resolution parity — a scale difference biases the peak location, so a
        #    shift measured on mismatched pixel sizes is not interpretable.
        if (abs(chm.res[0] - expected_res_m) > 1e-6
                or abs(chm.res[0] - sar.res[0]) > 1e-6):
            log.error("resolution_mismatch", chm_res=chm.res, sar_res=sar.res)
            raise ValueError("CHM and SAR must be resampled to a common resolution")

        # 3. Read the shared overlap window from both rasters.
        left = max(chm.bounds.left, sar.bounds.left)
        bottom = max(chm.bounds.bottom, sar.bounds.bottom)
        right = min(chm.bounds.right, sar.bounds.right)
        top = min(chm.bounds.top, sar.bounds.top)
        if right <= left or top <= bottom:
            raise ValueError("No spatial overlap between CHM and SAR")

        chm_win = chm.read(1, window=from_bounds(left, bottom, right, top, chm.transform))
        sar_win = sar.read(1, window=from_bounds(left, bottom, right, top, sar.transform))

    # 4. Trim to a common shape and standardise; NaNs corrupt the FFT peak.
    rows = min(chm_win.shape[0], sar_win.shape[0])
    cols = min(chm_win.shape[1], sar_win.shape[1])
    ref = np.nan_to_num(chm_win[:rows, :cols].astype("float64"))
    mov = np.nan_to_num(sar_win[:rows, :cols].astype("float64"))
    ref = (ref - ref.mean()) / (ref.std() + 1e-9)
    mov = (mov - mov.mean()) / (mov.std() + 1e-9)

    # 5. Subpixel phase cross-correlation: (dy, dx) aligns `mov` onto `ref`.
    shift, error, _ = phase_cross_correlation(ref, mov, upsample_factor=UPSAMPLE_FACTOR)
    dy, dx = float(shift[0]), float(shift[1])
    magnitude_px = float(np.hypot(dy, dx))
    drift_flagged = magnitude_px > DRIFT_TOLERANCE_PX

    report = {
        "shift_dy_px": round(dy, 3),
        "shift_dx_px": round(dx, 3),
        "shift_dy_m": round(dy * expected_res_m, 3),
        "shift_dx_m": round(dx * expected_res_m, 3),
        "magnitude_px": round(magnitude_px, 3),
        "phase_error": round(float(error), 4),
        "tolerance_px": DRIFT_TOLERANCE_PX,
        "drift_flagged": drift_flagged,
    }
    log.info("coregistration_drift_estimated", **report)
    return report

Two guards keep this diagnostic honest. A high phase_error — the normalised RMS of the correlation peak — signals that the two layers share little structure in this window (deep shadow, water, or genuine terrain distortion), so the shift estimate should not be trusted even if its magnitude looks small. And because the estimator can only recover a coherent translation, a tile that reports low drift here but still fuses poorly downstream is a strong indicator of the spatially variable terrain distortion described above, which must be routed back to geocoding rather than corrected with a shift.

Deterministic Transformation Logic #

Once drift is measured and attributed to a coherent offset, the correction is deterministic: reproject both layers onto a common equal-area grid so that later area-weighted biomass sums are honest, apply the estimated shift to the SAR layer as an affine translation, and re-verify that alignment actually improved. Using an equal-area CRS such as EPSG:6933 rather than the native UTM keeps per-pixel area constant across the tile, which matters because the fused product is ultimately integrated to tonnes. The correction is gated: it must both raise the normalised cross-correlation and clear an absolute floor, otherwise the tile is rejected rather than silently shipped with a shift that did nothing.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.transform import Affine
import structlog

log = structlog.get_logger()

EQUAL_AREA_CRS = "EPSG:6933"   # World Cylindrical Equal Area — area-true for AGB sums
CORR_GATE = 0.60               # minimum normalised cross-correlation after correction


def _reproject_to_grid(src_path: str, dst_crs: str, res_m: float):
    """Reproject band 1 of a raster onto a fresh grid at the target resolution."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=res_m)
        dst = np.empty((height, width), dtype="float32")
        reproject(source=rasterio.band(src, 1), destination=dst,
                  src_transform=src.transform, src_crs=src.crs,
                  dst_transform=transform, dst_crs=dst_crs,
                  resampling=Resampling.cubic)
    return dst, transform


def _ncc(a: np.ndarray, b: np.ndarray) -> float:
    """Normalised cross-correlation over the common footprint."""
    a = np.nan_to_num(a.astype("float64"))
    b = np.nan_to_num(b.astype("float64"))
    a = (a - a.mean()) / (a.std() + 1e-9)
    b = (b - b.mean()) / (b.std() + 1e-9)
    return float((a * b).mean())


def correct_coregistration(
    chm_path: str, sar_path: str, shift_dy_dx: tuple[float, float],
    out_path: str, dst_crs: str = EQUAL_AREA_CRS, res_m: float = 10.0,
) -> dict:
    """Reproject both layers to a common equal-area grid, apply the estimated
    affine translation to SAR, and gate on an improved cross-correlation."""
    chm, chm_tf = _reproject_to_grid(chm_path, dst_crs, res_m)
    sar, sar_tf = _reproject_to_grid(sar_path, dst_crs, res_m)
    log.info("reprojected_to_common_grid", crs=dst_crs, res_m=res_m,
             chm_shape=chm.shape, sar_shape=sar.shape)

    rows, cols = min(chm.shape[0], sar.shape[0]), min(chm.shape[1], sar.shape[1])
    chm, sar = chm[:rows, :cols], sar[:rows, :cols]
    ncc_before = _ncc(chm, sar)

    # Apply the (dy, dx) shift as an affine translation on the SAR grid, then
    # resample back onto the CHM grid so both bands share one transform.
    dy, dx = shift_dy_dx
    corrected_tf = sar_tf * Affine.translation(dx, dy)
    shifted = np.empty_like(sar)
    reproject(source=sar, destination=shifted,
              src_transform=corrected_tf, src_crs=dst_crs,
              dst_transform=chm_tf, dst_crs=dst_crs,
              resampling=Resampling.cubic)
    ncc_after = _ncc(chm, shifted)

    # Gate: the correction must materially improve alignment AND clear the floor.
    if ncc_after < ncc_before or ncc_after < CORR_GATE:
        log.error("coregistration_gate_failed", ncc_before=round(ncc_before, 4),
                  ncc_after=round(ncc_after, 4), gate=CORR_GATE)
        raise RuntimeError("Correction did not improve alignment past the gate")

    profile = {"driver": "GTiff", "height": rows, "width": cols, "count": 2,
               "dtype": "float32", "crs": dst_crs, "transform": chm_tf}
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(chm.astype("float32"), 1)     # band 1: reference CHM
        dst.write(shifted.astype("float32"), 2)  # band 2: corrected SAR
    log.info("coregistration_corrected", out_path=out_path,
             shift_dy=dy, shift_dx=dx,
             ncc_before=round(ncc_before, 4), ncc_after=round(ncc_after, 4))
    return {"ncc_before": ncc_before, "ncc_after": ncc_after, "out_path": out_path}

The Affine.translation(dx, dy) step is where discipline pays off: applying the shift to the SAR transform and then resampling onto the CHM grid guarantees both bands end up on one identical geotransform, so the fused stack is self-consistent by construction. Cubic resampling preserves the backscatter gradients the regression reads; nearest-neighbour would reintroduce edge aliasing that mimics the very drift being corrected.

Compliance Gating & Audit Trail Generation #

A correction that is not recorded is indistinguishable from tampering under third-party review. Every coregistration run therefore emits an append-only record binding the measured shift, the before-and-after cross-correlation, and the resampling parameters to the input checksums, so a verifier can reconstruct exactly what was moved and by how much. Under ISO 14064-3 this is the geometric-consistency evidence that underpins any positional-accuracy claim; under Verra VM-series methodologies it substantiates that the AGB surface was built on spatially reconciled inputs; and where the shift is material it must be traceable so that its downstream effect on reported tonnage is auditable rather than hidden.

The gate must fail loudly. A tile whose corrected cross-correlation clears the floor is written with its transform stamped into metadata and forwarded; a tile that cannot be aligned by translation is rejected and escalated to terrain re-correction, never coerced. That escalation flag, and the residual shift on accepted tiles, both feed the uncertainty budget consumed by emission factor uncertainty mapping, so that residual misregistration is propagated as variance rather than pretended away.

Production Integration #

In production the corrective loop runs as a fixed, containerised sequence with pinned dependencies, so every coregistration decision is reproducible byte-for-byte. Each tile follows the same six-step contract:

1. Ingest — pull the CHM and SAR products by identifier, verify SHA-256 checksums, and record the native CRS, grid origin, and resolution of each layer as source provenance.
2. Diagnose — enforce CRS and resolution parity, read the shared overlap window, and run estimate_coregistration_drift; admit tiles within tolerance directly and flag the rest for correction.
3. Transform — for flagged tiles, run correct_coregistration to reproject both layers onto the common equal-area grid and apply the estimated affine translation to the SAR band.
4. Validate — re-estimate cross-correlation after the shift and gate on both a material improvement and the absolute floor; reject and escalate any tile that fails, routing suspected terrain distortion back to range-Doppler geocoding.
5. Export — write the coregistered two-band stack to versioned, immutable object storage with the correction transform, tolerance, and gate metrics stamped into metadata.
6. Submit — attach the signed coregistration audit record and forward the validated stack to the fusion regression and the registry submission queue.

Batch execution chunks tiles through a scheduler, but the contract is per-tile so that a single drifted scene can never contaminate a neighbour. Read SAR windowed to the LiDAR footprint rather than whole-scene to keep memory bounded across thousands of tiles, and cache the measured shift alongside the tile identifier so re-runs skip re-estimation when inputs are unchanged. With this loop in place, coregistration drift stops being an invisible tax on biomass RMSE and becomes a measured, corrected, and audited property of every fused surface.

Related guides #

• Biomass Estimation from LiDAR–SAR Fusion — the parent stage whose fused stack this correction feeds.
• Ground-Truth Alignment for Carbon Models — plot-to-pixel calibration that requires a coregistered stack first.
• Debugging Silent Datum Shifts in Carbon Pipelines — distinguishing a transform bug from physical drift.
• Emission Factor Uncertainty Mapping — where residual misregistration is propagated as variance.