Harmonising Landsat and Sentinel-2 Reflectance

Pair bands by wavelength, put both on one grid, and apply a per-band linear adjustment to one sensor:

PAIRS = {"blue": ("SR_B2", "B02"), "green": ("SR_B3", "B03"), "red": ("SR_B4", "B04"),
         "nir": ("SR_B5", "B8A"), "swir1": ("SR_B6", "B11"), "swir2": ("SR_B7", "B12")}

harmonised = {k: landsat[l] * slope[k] + intercept[k] for k, (l, _) in PAIRS.items()}

Band 8A rather than 8 for the near infrared is the detail that removes most of the residual. This page belongs to atmospheric correction and surface reflectance in Satellite Processing Workflows & Index Pipelines.

Why the Bands Need Pairing Carefully

Band positions along the spectrum Visible bands of the two sensors overlap closely. In the near infrared, Sentinel-2 band 8 is much broader than Landsat's band 5, while Sentinel-2 band 8A is narrow and sits almost exactly where Landsat band 5 does. The shortwave infrared bands pair closely. Pairing on 8A rather than 8 is what makes the near-infrared harmonisation work. Near infrared is where pairing matters Landsat Sentinel-2 B8 broad B8A blue red NIR SWIR B8A sits almost exactly on Landsat's NIR; B8 spans far more of the spectrum.

The broad band 8 is the band most Sentinel-2 workflows use at 10 m, because 8A is 20 m. That is fine within Sentinel-2, but pairing it with Landsat’s narrow near-infrared band leaves a surface-dependent offset that no single linear coefficient fully removes. For harmonised 30 m products, use 8A.

Environment & Setup

Package Version pin Used for
rioxarray >=0.15 Loading and aligning both sensors
numpy >=1.23 Applying coefficients
scikit-learn >=1.3 Fitting local coefficients robustly 
pip install "rioxarray>=0.15" "numpy>=1.23" "scikit-learn>=1.3"

Complete Working Example

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

PAIRS = {"blue": ("SR_B2", "B02"), "green": ("SR_B3", "B03"), "red": ("SR_B4", "B04"),
         "nir": ("SR_B5", "B8A"), "swir1": ("SR_B6", "B11"), "swir2": ("SR_B7", "B12")}


def fit_coefficients(landsat: xr.Dataset, s2: xr.Dataset,
                     max_samples: int = 200_000, seed: int = 0) -> dict[str, tuple[float, float]]:
    """Per-band slope and intercept mapping Landsat onto Sentinel-2, from a same-day pair."""
    rng = np.random.default_rng(seed)
    coeffs = {}
    for name, (l_band, s_band) in PAIRS.items():
        x, y = landsat[l_band].values.ravel(), s2[s_band].values.ravel()
        ok = np.isfinite(x) & np.isfinite(y) & (x > 0) & (y > 0)
        idx = np.flatnonzero(ok)
        idx = rng.choice(idx, min(max_samples, idx.size), replace=False)
        m = HuberRegressor().fit(x[idx, None], y[idx])
        coeffs[name] = (float(m.coef_[0]), float(m.intercept_))
    return coeffs


def harmonise_landsat(landsat: xr.Dataset, coeffs: dict[str, tuple[float, float]]) -> xr.Dataset:
    out = {}
    for name, (l_band, _) in PAIRS.items():
        slope, intercept = coeffs[name]
        out[name] = (landsat[l_band] * slope + intercept).astype("float32")
    ds = xr.Dataset(out)
    ds.attrs.update(harmonised_to="sentinel-2", sensor="landsat")
    return ds

Fitting on a same-day pair — two acquisitions within a day of each other over the same ground — is what isolates the sensor difference from surface change. Both inputs must already be on one grid; the alignment is covered in matching Landsat and Sentinel-2 grids. The robust regressor matters because residual cloud and shadow produce asymmetric outliers that pull an ordinary least-squares line.

Published or Local Coefficients

Global coefficients or local ones Published coefficients from harmonisation studies are fitted over many sites worldwide and work well on average. Locally fitted coefficients from a few same-day pairs in your study area adapt to its surfaces and atmosphere and usually reduce residual differences further, at the cost of needing clear coincident acquisitions. Two sources of coefficients published, global fitted across many sites no coincident pairs needed good on average the default starting point fitted locally from same-day pairs in your area adapts to local surfaces usually smaller residuals worth it for unusual landscapes Fit locally, then compare with the published values — large disagreement is a warning sign.

Published coefficients are a sound default and require nothing but the numbers. Local coefficients are worth fitting when the study area’s surfaces are unlike the global average — arid, snow-covered, or dominated by one crop — and when a few clear same-day pairs are available. Comparing the two is informative in itself: close agreement confirms the fit, while a large disagreement usually means the local pair was contaminated by cloud or surface change.

After Harmonisation

Harmonised reflectance lets the two sensors share one time series, which roughly doubles the observation density — valuable wherever cloud limits the usable dates. It does not make the sensors identical. Landsat’s 30 m pixels are mixtures of what Sentinel-2’s 10 m pixels see separately, so harmonised products are best built at 30 m; view angles differ between overpasses; and the acquisition times differ by tens of minutes, so fast-changing surfaces such as flood water can genuinely differ between a same-day pair.

Keeping a sensor identifier on every observation preserves the ability to check any of those effects later, and it lets a model or a composite treat the sensors differently if the analysis requires it. The model-side view of the same problem is covered in harmonising features across sensors for transfer.

Verification

Residuals on a held-out pair On a same-day pair not used for fitting, the median absolute difference between Landsat and Sentinel-2 falls in every band after harmonisation, most in the near infrared where the band difference was largest. Residuals that do not fall indicate a band pairing or alignment problem. Median |Landsat − S2| on a held-out pair red green NIR before after Validate on a pair the coefficients were not fitted from. 
import numpy as np

held_out = harmonise_landsat(landsat_pair2, coeffs)
for name, (_, s_band) in PAIRS.items():
    before = np.nanmedian(np.abs(landsat_pair2[PAIRS[name][0]] - s2_pair2[s_band]))
    after = np.nanmedian(np.abs(held_out[name] - s2_pair2[s_band]))
    print(f"{name:<6} before {float(before):.4f}  after {float(after):.4f}")
    assert after <= before, f"harmonisation made {name} worse"

Validating on a pair the coefficients were not fitted from is the honest test, for the same reason a model is never evaluated on its training data. If residuals do not fall for a band, the pairing for that band is wrong or the two inputs are not properly aligned.

Common Errors

NIR residuals stay large after harmonisation

Sentinel-2 band 8 was paired with Landsat NIR. Use band 8A, resampled to the common 30 m grid alongside the other bands.

Coefficients vary wildly between pairs

The pairs include cloud, shadow or real surface change. Mask thoroughly and use pairs within a day of each other.

A sawtooth remains in the combined series

Coefficients were applied to the wrong sensor, or with slope and intercept swapped. Check that harmonised Landsat values move toward Sentinel-2 on a known pair.

Harmonised values fall outside 0–1

Coefficients fitted on scaled integers were applied to reflectance. Fit and apply on the same scale. Convert both sensors to physical reflectance first, then fit, then apply.

Frequently Asked Questions

Q: Which Sentinel-2 NIR band matches Landsat’s? Band 8A. Sentinel-2’s broad band 8 is much wider than Landsat’s near-infrared band and spans a water vapour feature, while the narrow band 8A is spectrally close to Landsat’s. Pairing band 8 with Landsat NIR leaves a systematic offset that coefficients only partly remove.

Q: Should Landsat be adjusted to Sentinel-2 or the reverse? Adjust whichever sensor contributes fewer observations to your series toward the one that contributes more, so most of the data is untouched. For recent periods that usually means adjusting Landsat to Sentinel-2.

Q: Is harmonisation enough to treat the sensors as one? For reflectance and indices at 30 metres, largely yes. Differences in spatial resolution, viewing geometry and revisit timing remain, so resample to a common resolution and record which sensor each observation came from.

Q: Are pre-harmonised products available? Harmonised Landsat and Sentinel-2 products exist and apply exactly these corrections, plus a common grid, at 30 metres. Where they cover your area and period, they save the work described here and are the natural choice for long, dense time series.

Q: Do the coefficients need refitting over time? Occasionally. Sensor calibration drifts slowly and providers issue collection updates, so coefficients fitted on pairs from one period may not suit another years later. Refit when a new collection or processing baseline is adopted, and keep the date range of the fitting pairs alongside the coefficients so their provenance is clear.