Harmonising Features across Sensors for Transfer

Map bands by wavelength, adjust with a linear fit, and check the result on overlapping acquisitions:

import numpy as np

# Per-band slope and intercept that put Landsat 8 onto the Sentinel-2 scale
COEFF = {"red": (0.982, 0.0012), "nir": (1.021, -0.0031), "swir1": (0.996, 0.0008)}

def harmonise(band: np.ndarray, name: str) -> np.ndarray:
    slope, intercept = COEFF[name]
    return (band * slope + intercept).astype("float32")

Without this step a model trained on one sensor learns to detect the sensor. This page belongs to feature engineering for pixel-based models in Raster Machine Learning & Model Inference.

Equivalent Bands Are Not Identical Bands

Two sensors described as having a “near infrared band” can differ by tens of nanometres in centre wavelength and by a factor of two in width. The surface is the same; the integral over the bandpass is not.

Two near infrared bands over the same spectrum A vegetation reflectance spectrum rises sharply at the red edge and plateaus through the near infrared with a water vapour absorption dip. One sensor's wide band spans the dip and returns a lower integrated reflectance; the other sensor's narrow band sits beside it and returns a higher value. The same leaf gives two different numbers. Same surface, two bandpasses wide NIR — spans the absorption dip narrow NIR water vapour absorption 650 nm 900 nm Integrating the green curve over the blue band gives a lower number than over the orange one.

Typical differences between Landsat 8 and Sentinel-2 surface reflectance are a few per cent per band — small in absolute terms and large relative to the decision boundaries a fitted model uses. An NDVI threshold tuned on one sensor can sit on the wrong side of a class boundary on the other.

Environment & Setup

Package Version pin Used for
numpy >=1.23 Coefficient application and the regression
scikit-learn >=1.3 Robust linear fits for locally derived coefficients
xarray >=2023.12 Keeping band names attached through the adjustment
rioxarray >=0.15 Loading matched pairs onto a common grid 
pip install "numpy>=1.23" "scikit-learn>=1.3" "xarray>=2023.12" "rioxarray>=0.15"

Complete Working Example

from dataclasses import dataclass

import numpy as np
import xarray as xr
from sklearn.linear_model import HuberRegressor


BAND_MAP = {          # secondary sensor band -> reference sensor band
    "SR_B4": "B04",   # red
    "SR_B5": "B08",   # near infrared
    "SR_B6": "B11",   # shortwave infrared 1
    "SR_B7": "B12",   # shortwave infrared 2
}


@dataclass
class Harmoniser:
    coefficients: dict[str, tuple[float, float]]

    def apply(self, ds: xr.Dataset) -> xr.Dataset:
        out = {}
        for src_name, ref_name in BAND_MAP.items():
            if src_name not in ds:
                continue
            slope, intercept = self.coefficients[ref_name]
            out[ref_name] = (ds[src_name] * slope + intercept).astype("float32")
        return xr.Dataset(out)

    @classmethod
    def fit(cls, secondary: xr.Dataset, reference: xr.Dataset,
            *, max_samples: int = 200_000, seed: int = 0) -> "Harmoniser":
        """Fit per-band coefficients from a near-simultaneous pair."""
        rng = np.random.default_rng(seed)
        coeffs: dict[str, tuple[float, float]] = {}

        for src_name, ref_name in BAND_MAP.items():
            x = secondary[src_name].values.ravel()
            y = reference[ref_name].values.ravel()
            ok = np.isfinite(x) & np.isfinite(y) & (x > 0) & (y > 0)
            idx = np.flatnonzero(ok)
            if idx.size > max_samples:
                idx = rng.choice(idx, max_samples, replace=False)

            # Huber rather than least squares: residual cloud is an outlier problem
            model = HuberRegressor(epsilon=1.35).fit(x[idx, None], y[idx])
            coeffs[ref_name] = (float(model.coef_[0]), float(model.intercept_))

        return cls(coeffs)

Two decisions in that fit matter. A robust regressor rather than ordinary least squares, because residual cloud and shadow are asymmetric outliers that drag a least-squares line noticeably. And fitting on a pair — two acquisitions within a day or two over the same ground — because anything longer apart mixes real surface change into the sensor difference.

Where Harmonisation Is Not Enough

What a linear adjustment can and cannot fix Gain and offset differences between sensors are exactly what a linear fit removes. Bandpass differences are partly removed, with a residual that depends on the surface type. Differences in spatial resolution and in acquisition geometry are not addressed at all by a reflectance adjustment and need resampling or a modelling decision instead. Three kinds of difference, one tool gain and offset calibration and atmospheric processing fully fixed by a linear fit bandpass shape different centre and width per band partly fixed; residual varies by surface resolution and geometry 30 m versus 10 m, different view angles not addressed at all A 30 m pixel is a mixture of what three 10 m pixels see — no reflectance coefficient recovers that. Either degrade the finer sensor to match, or train on the coarser one and accept the loss.

The resolution column is the one that surprises people. A model trained at 10 m on Sentinel-2 and applied to 30 m Landsat sees pixels that are spectral mixtures rather than pure surfaces, so its decision boundaries — tuned on relatively pure pixels — sit in the wrong place regardless of how well the reflectance was harmonised. The honest options are to train at the coarser resolution, or to keep separate models. Resampling Landsat up to 10 m does not create information, though it does at least make the grids match; the mechanics are in matching Landsat and Sentinel-2 grids.

A pragmatic middle path is to build the feature stack entirely from indices and temporal statistics rather than raw reflectance. Indices cancel much of the gain difference, and seasonal percentiles are robust to the residual, so a model built that way transfers far better than one keyed to absolute band values — at the cost of discarding whatever absolute brightness was telling you.

Verification

Transfer is a claim that has to be tested, and the test is not the training score.

import numpy as np

# 1. Distribution overlap: harmonised secondary should match the reference
for band in ("B04", "B08", "B11"):
    a = np.nanpercentile(reference[band], [10, 50, 90])
    b = np.nanpercentile(harmonised[band], [10, 50, 90])
    print(f"{band}: reference {a.round(4)}  harmonised {b.round(4)}")
    assert np.allclose(a, b, atol=0.02), f"{band} still offset after harmonisation"

# 2. Decision agreement on the overlap
agree = (model.predict(features_ref) == model.predict(features_harmonised)).mean()
print(f"class agreement across sensors: {agree:.1%}")
Before and after the linear adjustment Plotting matched pixels from two sensors against each other, the raw pairs sit on a line clearly offset from the one-to-one diagonal. After the fitted slope and intercept are applied the cloud of points straddles the diagonal, and the remaining spread is the irreducible bandpass and geometry difference rather than a systematic bias. raw pairs after harmonisation systematically above the diagonal straddling the diagonal; spread is irreducible

Anything below about 90% class agreement on a same-day overlap means the transfer is not working, and the cause is almost always either a band mapped to the wrong wavelength or a resolution mismatch that no coefficient can fix. Run this check on a held-out pair rather than on the pair the coefficients were fitted from, for the same reason any other model is not evaluated on its training data.

Common Errors

Coefficients near 1.0 but predictions still differ wildly

The bands were mapped by index rather than by wavelength — Landsat’s SR_B5 is near infrared while Sentinel-2’s band 5 is a red-edge band. Map by spectral position, always.

Harmonised values drift outside 0 to 1

The intercept was fitted on scaled integers and applied to reflectance, or vice versa. Fit and apply on the same scale, and assert the output range afterwards.

Transfer works in one region and fails in another

Globally published coefficients were fitted over different surfaces from yours. Refit locally on your own overlapping pairs, which is usually a few hours of work and worth it.

Frequently Asked Questions

Q: Why can a model not just be applied to another sensor directly? Because equivalent bands are not identical. Landsat 8’s near infrared is narrower than Sentinel-2’s and avoids a water vapour absorption feature, so the same surface returns a measurably different reflectance. A model with a threshold learned on one sensor applies that threshold to a shifted distribution on the other.

Q: Are indices more transferable than raw bands? Generally yes, because a ratio cancels multiplicative differences in gain and much of the illumination effect. It does not cancel bandpass differences, so an index still shifts between sensors — just less than the bands it is built from.

Q: How do I get harmonisation coefficients for my own data? Fit them on near-simultaneous acquisitions over your own study area: pair pixels from both sensors within a day or two, filter to cloud-free stable targets, and regress one on the other per band. Locally fitted coefficients beat published global ones when the landscape is unusual.