Spatial Pipeline EngineeringMRV Automation

Harmonizing Sentinel-2 and Landsat Surface Reflectance

Combining Sentinel-2 and Landsat into a single time series roughly doubles the observation density over a project area, which is the difference between a usable change signal and a cloud-shredded one in most tropical settings. It also introduces a systematic step in the data every time the series switches sensor, and that step has repeatedly been reported as forest degradation by pipelines that did not correct for it. This guide covers the corrections that make the combined series safe to difference, within Sentinel-2 and Landsat cloud masking workflows in the satellite imagery processing for emissions tracking stack.

Both sensors deliver a product labelled surface reflectance, and the label conceals four separate differences: the wavelengths each band actually integrates, the viewing geometry at acquisition, the ground sample distance, and the atmospheric correction applied. Each contributes an offset of a few percent in reflectance units — small individually, and collectively large enough to swamp the signal a degradation monitor is looking for, since canopy thinning also moves reflectance by a few percent.

A combined time series before and after harmonisation Two stacked plots of a vegetation index over three years. The upper plot shows the raw combined series: Sentinel-2 observations sit consistently above Landsat observations by a small but visible offset, so the series appears to step up and down as the sensor alternates, and a change detection algorithm fits breakpoints at the sensor switches. The lower plot shows the same observations after band pass adjustment and BRDF normalisation: the two sensors now overlap within their noise, the series is smooth, and the only breakpoint left is a real disturbance in the second year. An annotation notes that the apparent breakpoints in the upper plot are entirely artefacts of sensor alternation, and that they cluster at whatever cadence the two revisit cycles happen to produce. The same observations, before and after cross-sensor correction high low high low Raw combined saw-tooth from sensor alternation; breakpoints fitted at every switch Harmonised sensors overlap within noise; one real break remains the only real disturbance

Root Cause Analysis

Four independent mechanisms produce the offset, and they need separating because only two of them are correctable by a simple linear adjustment.

Band pass mismatch. Sentinel-2’s red band spans a different wavelength interval from Landsat 8’s, and its near-infrared band is markedly narrower. Because vegetation reflectance varies steeply across those intervals, the two instruments integrate different amounts of signal over what is nominally the same band. The resulting difference is systematic, depends on the surface being observed, and is largest exactly where the spectral response is steepest — the red edge, which is where vegetation monitoring lives. Published band pass adjustment coefficients handle the first-order part of this well; they are surface-type dependent, and applying coefficients derived over cropland to closed tropical forest leaves a residual.

Viewing and illumination geometry. Landsat views near-nadir; Sentinel-2 has a field of view of about twenty degrees, so an observation near the swath edge sees the canopy at a substantially different angle from one at the centre. Forest canopies are strongly anisotropic reflectors — they look brighter when viewed from the illumination direction and darker from the opposite side — so the same stand yields different reflectance depending on where in the swath it falls. This is a within-Sentinel-2 problem as well as a cross-sensor one, and it is the component most often left uncorrected.

Spatial support. Ten and twenty metre Sentinel-2 pixels against thirty metre Landsat pixels means each observation integrates a different mix of canopy, gap, and shadow. Resampling changes the numbers but not the underlying support mismatch, and in heterogeneous canopy the residual difference persists after every other correction. This one cannot be removed, only managed, by aggregating both to a common coarser grid where support genuinely matches.

Atmospheric correction. The two agencies’ processors differ in their aerosol retrieval and their handling of adjacency effects, and the difference is largest under high aerosol loading — which in practice means the tropics in burning season, precisely when a deforestation monitor most wants data. Product version matters here too: reprocessing campaigns have shifted these offsets, so a series spanning a collection change contains a step even within a single sensor.

The pattern to internalise is that the first two are correctable with published coefficients and a BRDF model, the third is manageable by choice of grid, and the fourth is best handled by filtering rather than correcting — excluding observations under high aerosol loading rather than trying to fix them.

Diagnostic Pipeline / Pre-Flight Validation

Before harmonising anything, measure the offset empirically on the project’s own surfaces. Published coefficients are a starting point, not an answer, and the check that matters is whether same-day observations from the two sensors agree after correction.

from dataclasses import dataclass
from datetime import date

import numpy as np
import structlog

log = structlog.get_logger()


@dataclass(frozen=True)
class CrossSensorPair:
    """Same location, same day, both sensors — the only clean comparison.

    Same-day pairs are rare (the orbits coincide occasionally) but they are
    the only observations where the surface is guaranteed identical, so
    every offset estimate should be built from them rather than from
    seasonal means that confound phenology with sensor difference.
    """
    lon: float
    lat: float
    acquired: date
    band: str
    s2_reflectance: float
    landsat_reflectance: float
    s2_view_zenith_deg: float
    solar_zenith_deg: float
    land_cover: str


@dataclass(frozen=True)
class BandOffset:
    band: str
    slope: float
    intercept: float
    residual_std: float
    n_pairs: int
    land_cover: str


MIN_PAIRS = 200
MAX_ACCEPTABLE_RESIDUAL = 0.015     # reflectance units


def fit_band_offset(pairs: list[CrossSensorPair], band: str, cover: str) -> BandOffset:
    """Ordinary least squares of Landsat on Sentinel-2 for one band and cover.

    Fitted per land cover deliberately. A single global adjustment fitted
    across forest, water, and bare soil is dominated by whichever class has
    the most pairs, and applies that class's correction everywhere.
    """
    sel = [p for p in pairs if p.band == band and p.land_cover == cover]
    if len(sel) < MIN_PAIRS:
        raise ValueError(
            f"{len(sel)} same-day pairs for band {band} over {cover}; "
            f"at least {MIN_PAIRS} are needed for a stable fit. Widen the "
            "date window to ±1 day before falling back to published "
            "coefficients, and record which was used."
        )

    x = np.array([p.s2_reflectance for p in sel])
    y = np.array([p.landsat_reflectance for p in sel])
    slope, intercept = np.polyfit(x, y, 1)
    residual = y - (slope * x + intercept)

    offset = BandOffset(
        band=band,
        slope=float(slope),
        intercept=float(intercept),
        residual_std=float(residual.std()),
        n_pairs=len(sel),
        land_cover=cover,
    )

    if offset.residual_std > MAX_ACCEPTABLE_RESIDUAL:
        log.warning(
            "harmonisation.residual_high",
            band=band, land_cover=cover,
            residual_std=round(offset.residual_std, 4),
            hint="check for view-angle spread — fit BRDF before band pass",
        )

    log.info(
        "harmonisation.offset_fitted",
        band=band, land_cover=cover,
        slope=round(offset.slope, 4), intercept=round(offset.intercept, 5),
        n_pairs=offset.n_pairs,
    )
    return offset


def assert_view_angle_spread(pairs: list[CrossSensorPair]) -> None:
    """Refuse to fit a band pass offset while view angle is uncorrected.

    If the Sentinel-2 pairs span a wide range of view zenith angles, the
    fitted 'band pass' offset silently absorbs the BRDF effect, and the
    resulting coefficients only work at the mean view angle of the fit set.
    """
    angles = np.array([p.s2_view_zenith_deg for p in pairs])
    spread = float(angles.max() - angles.min())
    if spread > 8.0:
        raise ValueError(
            f"view zenith spread is {spread:.1f}° across the fit set — apply "
            "BRDF normalisation to nadir first, then fit the band pass "
            "offset on normalised reflectance"
        )

The view-angle assertion is the one that catches the subtle version of this problem. A band pass correction fitted on an angle-diverse set looks excellent on its own fit data and then performs poorly on new observations acquired at a different position in the swath, which shows up months later as seasonal-looking drift.

The correction order, and what goes wrong when two steps are swapped A pipeline in five ordered steps with a warning attached to one transition. First, mask cloud and cloud shadow, because contaminated pixels corrupt every subsequent fit. Second, filter on aerosol optical depth, discarding rather than correcting observations under heavy haze. Third, normalise for bidirectional reflectance to a common nadir view and fixed solar zenith, which removes the angular component. Fourth, apply band pass adjustment per land cover on the normalised reflectance. Fifth, aggregate both sensors to a common grid so spatial support matches. A red annotation between steps three and four notes that swapping them causes the band pass coefficients to absorb the angular effect, producing coefficients valid only at the mean view angle of the fitting set. Order matters — steps 3 and 4 are not commutative 1 · Mask cloud and shadow contamination breaks every fit 2 · Filter on aerosol depth discard, do not correct 3 · BRDF to nadir fixed solar zenith, removes the angular term 4 · Band pass per land cover, on normalised reflectance 5 · Common grid match spatial support, not just resolution Swap 3 and 4 and the band pass coefficients absorb the angular effect. They then fit their own training set beautifully and fail on any observation acquired elsewhere in the swath — which surfaces later as a seasonal-looking drift with no physical explanation. the two amber steps must run in this order

Deterministic Transformation Logic

The harmonisation itself applies the corrections in the order established above and records what it did to each observation. The record matters as much as the correction, because a downstream analyst seeing an unexpected value needs to know which coefficients touched it.

from dataclasses import dataclass, replace


@dataclass(frozen=True)
class Observation:
    obs_id: str
    sensor: str                 # S2A | S2B | L8 | L9
    band: str
    reflectance: float
    view_zenith_deg: float
    solar_zenith_deg: float
    relative_azimuth_deg: float
    land_cover: str
    aerosol_optical_depth: float
    harmonised: bool = False
    corrections: tuple[str, ...] = ()


# c-factor BRDF normalisation, RossThick-LiSparse kernel weights per band.
# These are the MODIS-derived constants in standard use for this correction.
BRDF_WEIGHTS = {
    "red":  (0.0409, 0.0071),
    "nir":  (0.1789, 0.0574),
    "swir": (0.2580, 0.0605),
}
AOD_LIMIT = 0.4
TARGET_SOLAR_ZENITH = 45.0


def _kernel_ratio(band: str, view_zenith: float, solar_zenith: float) -> float:
    """Ratio of the kernel-weighted BRDF at nadir to that at the observation.

    A simplified c-factor. In production this calls the full RossThick and
    LiSparse kernels; the structure is what matters here — the correction is
    a multiplicative ratio, so it is exactly reversible and can be logged as
    a single number per observation.
    """
    import math

    f_iso = 1.0
    f_vol, f_geo = BRDF_WEIGHTS[band]
    vz, sz = math.radians(view_zenith), math.radians(solar_zenith)
    tz = math.radians(TARGET_SOLAR_ZENITH)

    at_obs = f_iso + f_vol * (vz + sz) * 0.35 + f_geo * (vz * sz) * 0.2
    at_nadir = f_iso + f_vol * tz * 0.35
    return at_nadir / at_obs


def harmonise(obs: Observation, offsets: dict[tuple[str, str], BandOffset]) -> Observation | None:
    """Bring one observation onto the Landsat reference scale.

    Returns None for observations that must be discarded rather than
    corrected. Landsat is chosen as the reference because its archive is
    longer, so harmonising toward it keeps historical series intact rather
    than requiring the whole archive to be reprocessed.
    """
    if obs.harmonised:
        raise ValueError(f"observation {obs.obs_id} has already been harmonised")

    if obs.aerosol_optical_depth > AOD_LIMIT:
        log.info("harmonisation.discarded", obs_id=obs.obs_id,
                 reason="aod_above_limit", aod=obs.aerosol_optical_depth)
        return None

    value = obs.reflectance
    applied: list[str] = []

    ratio = _kernel_ratio(obs.band, obs.view_zenith_deg, obs.solar_zenith_deg)
    value *= ratio
    applied.append(f"brdf_cfactor={ratio:.5f}")

    if obs.sensor.startswith("S2"):
        key = (obs.band, obs.land_cover)
        if key not in offsets:
            log.warning("harmonisation.no_offset", obs_id=obs.obs_id,
                        band=obs.band, land_cover=obs.land_cover)
            return None
        off = offsets[key]
        value = off.slope * value + off.intercept
        applied.append(f"bandpass={off.slope:.4f}x+{off.intercept:.5f}")

    if not 0.0 <= value <= 1.0:
        log.warning("harmonisation.out_of_range", obs_id=obs.obs_id,
                    corrected=round(value, 4), original=obs.reflectance)
        return None

    return replace(obs, reflectance=value, harmonised=True,
                   corrections=tuple(applied))

Two details in that function carry more weight than their line count. Refusing to harmonise an already-harmonised observation prevents the double-correction that occurs whenever a backfill overlaps a completed range, and it fails loudly instead of producing a plausible wrong number. Returning None rather than a flagged value for discarded observations means a caller cannot accidentally include them, which a nullable field invites.

Compliance Gating & Audit Trail Generation

A harmonised series carries an obligation that a single-sensor one does not: to show that any detected change is not a sensor artefact.

Per-observation correction records. The corrections tuple above is the audit unit — it states the exact multiplicative and linear factors applied, so any value can be reversed to its source.

The fitted offsets with their provenance. Which land covers were fitted, how many same-day pairs supported each, the residual standard deviation, and whether published coefficients were substituted for any band. A verifier comparing two monitoring periods will want to know whether the coefficients changed between them.

A sensor-composition summary per period. If one period is eighty percent Sentinel-2 and the next is forty percent, any residual harmonisation error appears as a change between the periods. Reporting the composition makes the risk visible and lets a reviewer weigh it.

Discard counts by reason. High aerosol, missing offset, out-of-range after correction — these should be counted and reported, because a period where a third of the observations were discarded for aerosol is a period with a thin, potentially biased sample, and the bias is toward clear days.

Production Integration

In practice the harmonisation sits between cloud masking and any temporal compositing, and it must run before change detection rather than after. The clearest sign it has been placed wrongly is breakpoints clustering at intervals matching the sensors’ combined revisit pattern, which the algorithms described in CCDC vs LandTrendr vs BFAST for carbon monitoring will happily fit as real events.

Where a project already consumes a harmonised product — HLS being the widely used one — most of this work is done upstream, and the remaining job is to verify rather than repeat it. The verification is the same same-day pair analysis: pull the pairs, fit the residual offset, and confirm it is within noise. Harmonised products are generally good, and they are fitted globally, so a project over an unusual surface can still see a residual worth knowing about.

The masking that precedes all of this deserves its own attention, since a cloud-contaminated pixel entering an offset fit corrupts the coefficients for every observation afterwards — see automating Sentinel-2 cloud masking with STAC and rasterio for that layer, and troubleshooting cloud shadow false positives in Sentinel-2 for the shadow cases that survive it.

Why the same canopy reflects differently across a Sentinel-2 swath A cross-section of a forest canopy with the sun at a fixed position on the left and three viewing positions across the swath. Viewing from the western edge, close to the illumination direction, the sensor sees mostly sunlit crown and little shadow, so reflectance reads high. Viewing at nadir in the swath centre, the sensor sees a mixture of sunlit crown and shadow, so reflectance reads at the reference level. Viewing from the eastern edge, away from the illumination direction, the sensor sees proportionally more shadow between crowns, so reflectance reads low. A panel notes that this variation across a single scene is comparable in size to the cross-sensor offset, and that it is the component most often left uncorrected. The same stand, three positions in the swath, three reflectances sun west edge mostly sunlit crown reflectance high nadir crown and shadow mixed reference level east edge more shadow between crowns reflectance low This within-scene variation rivals the cross-sensor offset — and it is the one most pipelines never correct.

Frequently Asked Questions

Should the series be harmonised toward Sentinel-2 or toward Landsat?

Toward Landsat, in almost every case, because the Landsat archive extends back four decades and a baseline period nearly always draws on it. Harmonising toward Sentinel-2 means every historical observation needs adjusting, which is both more computation and more risk, and it makes the series depend on a mission that started in 2015. The exception is a project whose entire history is post-2017 and whose analysis is resolution-sensitive, where keeping the ten-metre detail may be worth more than the archive.

Are published band pass coefficients good enough on their own?

For most land covers, close. They are fitted globally, so they carry a residual over any specific surface, and that residual is largest over surfaces poorly represented in the fitting set — dense tropical canopy, wetlands, and anything with strong understory contribution. The practical approach is to apply published coefficients as the default, fit local ones where enough same-day pairs exist, and record which was used per band and land cover. A project that never checks will not know which case it is in.

How many same-day pairs realistically exist?

Fewer than intuition suggests, and heavily biased toward clear conditions. Landsat and Sentinel-2 orbits coincide at a given location only occasionally, and both must be cloud-free at that moment. Over a single project area a year of data might yield a few hundred usable pixel pairs, which is why widening to a one-day window is common. The widening is safe for reflectance over stable surfaces and unsafe during rapid phenological change, so it should be applied outside the green-up and senescence windows.

Does harmonisation help or hurt a change detection algorithm’s sensitivity?

It helps substantially, and the mechanism is worth understanding. Change detection algorithms estimate a noise level from the series and flag departures beyond it. An unharmonised series has an inflated noise level because of the sensor saw-tooth, so genuine small changes fall inside the noise and go undetected. Harmonising both removes the false breakpoints at sensor switches and lowers the noise floor, making real degradation detectable that previously was not.

What about Landsat 7, and the scan line corrector gap?

Landsat 7 can be included, with two caveats. Its post-2003 imagery has systematic data gaps from the scan line corrector failure, which must be masked rather than interpolated for change work, and its radiometry differs enough from Landsat 8 that it needs its own offset rather than being pooled into a generic Landsat class. Many projects exclude it for the monitoring period and use it only for the historical baseline, where its lower quality is offset by there being no alternative.

How should a collection or processing baseline change be handled mid-series?

As a new sensor. A reprocessing campaign that changes the atmospheric correction produces a step in the series exactly like a cross-sensor one, and treating collection version as part of the sensor identity handles it with machinery already present. This means the offset table is keyed on sensor, collection, band, and land cover — larger, but it makes a collection migration a matter of fitting one more set of coefficients rather than discovering a mysterious step months afterward.

Is aggregating both sensors to a coarser common grid always necessary?

Not always, but it is the honest option when the analysis differences two dates. Resampling Sentinel-2 to thirty metres and calling the support matched is not quite true — a resampled ten-metre pixel and a native thirty-metre pixel integrate differently — though the residual is small relative to the other terms in closed canopy. In heterogeneous or fragmented landscapes the residual grows, and aggregating both to sixty metres or to a stand-level polygon removes it entirely at the cost of spatial detail. Which trade is right depends on whether the minimum mapping unit is closer to a pixel or a stand.