Matching Landsat and Sentinel-2 Grids

To combine the two missions, choose one target grid and resample both onto it — downsampling Sentinel-2 rather than upsampling Landsat, in most cases:

import rioxarray
from rasterio.enums import Resampling

landsat = rioxarray.open_rasterio("LC09_B5.tif", masked=True).squeeze(drop=True)      # 30 m
sentinel = rioxarray.open_rasterio("S2_B8A.tif", masked=True).squeeze(drop=True)      # 20 m

# Landsat's grid is the reference; average preserves radiometry when downsampling
sentinel_30m = sentinel.rio.reproject_match(landsat, resampling=Resampling.average)

This is the cross-mission case of the resampling decisions set out in Advanced Resampling and Upscaling Techniques.

Why This Arises in Remote Sensing Workflows

Combining Landsat and Sentinel-2 roughly triples the number of clear observations available for a site, which is transformative for time-series work in cloudy regions and for any application needing a short revisit. Harmonised products exist, but plenty of pipelines still have to do the combination themselves — because the harmonised product does not cover the area, does not include the band needed, or lags the raw archive by too long.

The mechanics are the same as any resampling problem, with three complications specific to this pair. The grids do not align, because the two missions use different tiling schemes even when both are in UTM. The resolutions differ by a factor that is not a clean integer for every band pairing. And the band centres differ, so even after perfect geometric alignment, the same index computed from each mission returns slightly different numbers over identical ground.

Only the first two are geometry. The third is radiometry, and no amount of resampling addresses it — it has to be measured and, if it matters, corrected.

Band correspondence between the two missions Blue, green, red and short-wave infrared bands correspond closely between Landsat 8/9 and Sentinel-2. The near-infrared is where they differ most: Sentinel-2 B08 is broad, B8A is narrow, and Landsat B5 sits between them, so an index using near-infrared is the one most affected by mission differences. Where the bands line up, and where they do not 0.49 µm 0.66 µm 0.86 µm 2.2 µm Sentinel-2 B02 B03 B04 B08 (wide) B8A B12 Landsat 8/9 B2 B3 B4 B5 B7 near-infrared: approximate, not equal Green lines are close correspondences; dashed lines are the pairing that needs a measured offset.

Environment & Setup

Package Version Why
rioxarray ≥0.15 reproject_match and masked reads
rasterio ≥1.3.0 Resampling enums, transforms
xarray ≥2023.6 Stacking the harmonised bands
numpy ≥1.23 Offset estimation 
pip install "rioxarray>=0.15" "rasterio>=1.3.0" "xarray>=2023.6"

Complete Working Example

This function builds a common-grid stack from one Landsat scene and one Sentinel-2 scene, estimates the cross-mission index offset over stable ground, and returns both the stack and the offset.

Native resolutions across the two missions Landsat delivers 15, 30 and resampled-100 metre bands while Sentinel-2 delivers 10, 20 and 60. Only the 30 metre optical set and the 20 metre Sentinel bands are close enough to combine without a resolution decision, which is why 30 metres is the usual common grid. What has to move to reach a common grid band group native resolution to a 30 m grid Landsat B2–B7 30 m unchanged — the reference Landsat B8 pan 15 m downsample with average Landsat B10/B11 100 m, delivered at 30 already there, but smooth Sentinel-2 B02–B04, B08 10 m downsample with average Sentinel-2 B8A, B11, B12 20 m downsample with average Every row but the first is a resampling step, and each one is an entry in the error budget. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling

# Pair by wavelength, not by band number
BAND_MAP = {
    "blue": {"landsat": "B2", "sentinel": "B02"},
    "green": {"landsat": "B3", "sentinel": "B03"},
    "red": {"landsat": "B4", "sentinel": "B04"},
    "nir": {"landsat": "B5", "sentinel": "B8A"},     # B8A, not B08: narrower, closer to L8 B5
    "swir22": {"landsat": "B7", "sentinel": "B12"},
}


def ndvi(nir: xr.DataArray, red: xr.DataArray) -> xr.DataArray:
    denom = nir + red
    return xr.where(denom != 0, (nir - red) / denom, np.nan).astype("float32")


def harmonise(
    landsat_paths: dict[str, str],
    sentinel_paths: dict[str, str],
    *,
    stable_mask: xr.DataArray | None = None,
) -> tuple[xr.Dataset, float]:
    """Resample Sentinel-2 onto the Landsat 30 m grid and measure the NDVI offset."""
    ls = {k: rioxarray.open_rasterio(p, masked=True).squeeze(drop=True).astype("float32")
          for k, p in landsat_paths.items()}
    reference = ls["red"]

    s2 = {}
    for key, path in sentinel_paths.items():
        da = rioxarray.open_rasterio(path, masked=True).squeeze(drop=True).astype("float32")
        # Downsampling reflectance: average, which preserves the mean over the footprint
        s2[key] = da.rio.reproject_match(reference, resampling=Resampling.average)

    ndvi_ls = ndvi(ls["nir"], ls["red"])
    ndvi_s2 = ndvi(s2["nir"], s2["red"])

    diff = ndvi_s2 - ndvi_ls
    if stable_mask is not None:
        diff = diff.where(stable_mask)
    offset = float(diff.median(skipna=True))

    stack = xr.Dataset({f"ls_{k}": v for k, v in ls.items()} |
                       {f"s2_{k}": v for k, v in s2.items()} |
                       {"ndvi_landsat": ndvi_ls, "ndvi_sentinel": ndvi_s2})
    stack.attrs["cross_mission_ndvi_offset"] = round(offset, 4)
    return stack, offset


if __name__ == "__main__":
    stack, offset = harmonise(
        {"red": "LC09_B4.tif", "nir": "LC09_B5.tif"},
        {"red": "S2_B04.tif", "nir": "S2_B8A.tif"},
    )
    print("median NDVI offset (S2 − Landsat):", offset)
    # Apply the offset if the two will be used in one series
    stack["ndvi_sentinel_adj"] = stack["ndvi_sentinel"] - offset

The offset is typically a few hundredths of an NDVI unit and is stable enough over a scene to be worth applying. Measuring it per scene pair rather than adopting a published constant is the more defensible route, because it absorbs atmospheric and illumination differences alongside the spectral ones.

Variant Patterns

1. When to upsample instead

Downsampling is the default, but a product that must align with an existing 10 m dataset has no choice. In that case upsample Landsat with bilinear and be explicit in the metadata that the effective resolution is still 30 m.

landsat_10m = landsat.rio.reproject_match(sentinel_10m_reference, resampling=Resampling.bilinear)
landsat_10m.attrs["effective_resolution_m"] = 30      # the honest statement

Any per-pixel statistic computed on the upsampled layer treats nine correlated pixels as independent observations, which inflates apparent precision — the reason the parent topic argues for resampling as few times as possible.

2. A shared grid neither mission owns

For multi-mission archives it is often cleanest to define the analysis grid once and match everything to it, so no mission is privileged and the grid does not change when a new sensor is added.

from rasterio.transform import from_origin

res = 30.0
grid = xr.DataArray(
    np.zeros((3000, 3000), dtype="float32"), dims=("y", "x"),
    coords={"y": 9_900_000 - (np.arange(3000) + 0.5) * res,
            "x": 300_000 + (np.arange(3000) + 0.5) * res},
)
grid.rio.write_crs("EPSG:32636", inplace=True)
grid.rio.write_transform(from_origin(300_000, 9_900_000, res, res), inplace=True)

Snapping the origin to a multiple of the resolution keeps tiles from different runs abutting exactly, as described in Aligning Two Rasters with reproject_match.

3. Quality layers across missions

Two quality encodings, one validity mask Landsat encodes quality as bit flags inside QA_PIXEL, while Sentinel-2 uses integer class codes in the SCL layer. Both must be translated into the same boolean convention before the missions can be combined, and both must be resampled with nearest. Normalise quality before combining anything else Landsat QA_PIXEL bit 3 cloud · bit 4 shadow bit 5 snow · bit 7 water Sentinel-2 SCL codes 3, 8, 9, 10, 11 invalid codes 4, 5, 6 keep boolean `clear` one convention, both missions resample nearest only Translate to boolean first, then resample: resampling bit flags or class codes with anything else corrupts both. 
import numpy as np

# Landsat: bits 3 (cloud) and 4 (cloud shadow) of QA_PIXEL
qa = rioxarray.open_rasterio("LC09_QA_PIXEL.tif").squeeze(drop=True).astype("uint16")
clear_ls = ((qa & (1 << 3)) == 0) & ((qa & (1 << 4)) == 0)

# Sentinel-2: SCL class codes
scl = rioxarray.open_rasterio("S2_SCL.tif").squeeze(drop=True)
clear_s2 = ~scl.isin([0, 1, 3, 8, 9, 10, 11])
clear_s2 = clear_s2.rio.reproject_match(clear_ls, resampling=Resampling.nearest)

clear_both = clear_ls & clear_s2

Verifying Co-registration

Geometric agreement between missions is the assumption most worth testing, because it fails quietly and its symptom — edge-aligned differences — looks like real change.

Difference the same index from both missions over a date pair a day or two apart, and inspect where the residual is largest. A well-registered pair shows a residual that is roughly uniform across the scene; a misregistered pair shows a residual concentrated along field boundaries, roads and coastlines, which is the signature described in Change Detection and Differencing Workflows.

If the residual is edge-aligned, the fix is a shift estimate rather than a different resampling method: cross-correlate a small window over a high-contrast area, find the offset that minimises the residual, and apply it before resampling. A half-pixel correction at 30 m is worth making; anything above one pixel points at a metadata problem rather than a registration one.

Common Errors

NDVI from the two missions differs by a constant

Expected: the band centres are not identical. Measure the offset over stable ground and apply it, rather than assuming the two are directly comparable.

The resampled Sentinel-2 layer looks blocky

nearest was used when downsampling reflectance, so each 30 m pixel takes a single 20 m sample instead of averaging the footprint. Use average.

Shapes match but the data is offset by one pixel

The two grids were assumed compatible rather than matched explicitly. Use reproject_match against one reference and assert transform equality.

Frequently Asked Questions

Q: Should I upsample Landsat to 10 m or downsample Sentinel-2 to 30 m? Downsample to 30 m for most analyses. Upsampling Landsat creates the appearance of detail that the sensor never captured, and any statistic computed on it double-counts information. Upsample only when the output must align with an existing 10 m product.

Q: Are Landsat 8 band 5 and Sentinel-2 band 8 interchangeable? They are both near-infrared but not identical: Sentinel-2 B08 is wide and B8A is narrow, and Landsat B5 sits between them. Indices computed from each differ by a small, systematic amount that should be measured over stable targets before the two are combined.

Q: Do the two missions share a grid origin? Not in general. Both are commonly delivered in UTM, but their tiling schemes and origins differ, so co-location requires an explicit resampling onto one chosen grid rather than an assumption that the pixels line up.