Resampling Sentinel-2 20m Bands to 10m

Sentinel-2 delivers B05, B11, and B12 at 20 m and B04/B08 at 10 m, so before you can compute an index across those bands they must share one grid. Open a native 10 m band as the template and reproject_match each 20 m band onto it with bilinear resampling:

import rioxarray
from rasterio.enums import Resampling

ref_10m = rioxarray.open_rasterio("B04.tif")     # native 10 m reference grid
swir_20m = rioxarray.open_rasterio("B11.tif")    # 20 m band to upsample

swir_10m = swir_20m.rio.reproject_match(
    ref_10m, resampling=Resampling.bilinear,
)
# swir_10m now shares B04's CRS, transform, width, and height exactly

This page sits within Handling Pixel Resolution and Scaling and covers the specific case of aligning Sentinel-2’s mixed 10 m / 20 m bands so downstream index math is pixel-for-pixel correct.

Why This Arises in Remote Sensing Workflows

Sentinel-2’s Multispectral Instrument samples different bands at different ground sampling distances: the visible and near-infrared bands (B02, B03, B04, B08) at 10 m; the red-edge and shortwave-infrared bands (B05, B06, B07, B8A, B11, B12) at 20 m; and atmospheric bands at 60 m. Any index that crosses those groups — NDMI from B08 and B11, a red-edge NDVI from B08 and B05, or a burn ratio from B08 and B12 — needs the operands on an identical grid.

If you simply load a 10 m band (shape 10980×10980) and a 20 m band (shape 5490×5490) and try to subtract them, the arrays do not broadcast and the operation fails. Force them to the same shape carelessly and you get a subtler bug: the pixel origins are offset by half a pixel, so band A’s pixel (0,0) does not cover the same ground as band B’s pixel (0,0). Per-pixel index values are then computed from mismatched footprints and are quietly wrong.

The correct fix is grid snapping: pick one native 10 m band as the reference and resample every coarser band to its exact transform. This is the alignment prerequisite that feeds spectral index calculation and any multi-band stacking. For the general theory of when to use each interpolation kernel, see Choosing the Right Resampling Method for Sentinel-2 in the Advanced Resampling & Upscaling Techniques cluster. For the broader context, see Core Raster Fundamentals & STAC Mapping.

Grid Snapping: 20 m Onto the 10 m Reference

The diagram shows why reproject_match matters: it aligns the coarse grid to the fine one so cell edges coincide, rather than floating on an offset origin.

Bilinear interpolation is the right kernel for continuous reflectance: each output 10 m value is a distance-weighted blend of the four nearest 20 m samples, which preserves the smooth radiometric surface. Nearest-neighbour would replicate whole 20 m blocks and reintroduce the coarse look; it belongs only on categorical layers such as the SCL mask, where blending class codes is meaningless.

The 60 m atmospheric bands can be upsampled the same way, but a 6× jump interpolates across a very coarse signal, so treat those results as approximate context rather than per-pixel truth.

Environment & Setup

pip install "rioxarray>=0.13.0" "rasterio>=1.3.0"

Quick check that the two grids differ before you align them:

import rioxarray

b04 = rioxarray.open_rasterio("B04.tif")   # 10 m
b11 = rioxarray.open_rasterio("B11.tif")   # 20 m
print("B04:", b04.rio.resolution(), b04.shape)
print("B11:", b11.rio.resolution(), b11.shape)  # coarser cells, smaller shape

Complete Working Example

This script aligns a set of 20 m bands to a 10 m reference, stacks them, and asserts pixel-perfect alignment before returning the cube.

import rioxarray
import xarray as xr
from rasterio.enums import Resampling


def align_s2_bands_to_10m(
    reference_10m_path: str,
    coarse_20m_paths: dict[str, str],
) -> xr.Dataset:
    """
    Upsample Sentinel-2 20 m bands onto a native 10 m reference grid.

    Parameters
    ----------
    reference_10m_path : str
        Path to a native 10 m band (B04 or B08) defining the target grid.
    coarse_20m_paths : dict[str, str]
        Mapping of band name -> path for each 20 m band to align.

    Returns
    -------
    xr.Dataset
        All bands on one 10 m grid, ready for index computation.
    """
    # The reference defines CRS, transform, width and height for everything.
    reference = rioxarray.open_rasterio(reference_10m_path, masked=True)

    aligned: dict[str, xr.DataArray] = {"reference": reference.squeeze()}

    for name, path in coarse_20m_paths.items():
        band_20m = rioxarray.open_rasterio(path, masked=True)
        # reproject_match copies the reference grid exactly onto the 20 m band.
        # bilinear = smooth interpolation appropriate for continuous reflectance.
        band_10m = band_20m.rio.reproject_match(
            reference, resampling=Resampling.bilinear,
        )
        aligned[name] = band_10m.squeeze()

    # Every array now shares one grid, so a Dataset merge is unambiguous.
    cube = xr.Dataset(aligned)

    # Alignment gate: shapes, transform and CRS must match the reference.
    ref_shape = reference.rio.shape
    ref_transform = reference.rio.transform()
    for name, arr in aligned.items():
        assert arr.rio.shape == ref_shape, f"{name} shape mismatch"
        assert arr.rio.transform() == ref_transform, f"{name} transform mismatch"
        assert arr.rio.crs == reference.rio.crs, f"{name} CRS mismatch"

    return cube


if __name__ == "__main__":
    cube = align_s2_bands_to_10m(
        reference_10m_path="B08.tif",           # 10 m NIR as the template
        coarse_20m_paths={"b05": "B05.tif",     # red-edge, 20 m
                          "b11": "B11.tif",     # SWIR-1, 20 m
                          "b12": "B12.tif"},    # SWIR-2, 20 m
    )
    # NDMI = (NIR - SWIR1) / (NIR + SWIR1), now valid pixel-for-pixel.
    nir, swir1 = cube["reference"], cube["b11"]
    ndmi = (nir - swir1) / (nir + swir1)
    print("NDMI grid:", ndmi.rio.shape, ndmi.rio.resolution())

Passing masked=True promotes the file’s nodata value to NaN, so bilinear interpolation does not smear a sentinel like 0 into valid pixels along scene edges.

Variant Patterns

1. Pure rasterio.warp — no xarray dependency

When you work in rasterio arrays rather than DataArrays, reproject gives the same grid snapping by passing the reference band’s transform, crs, and shape as the destination:

import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling

with rasterio.open("B04.tif") as ref:          # 10 m template
    dst_profile = ref.profile
    dst = np.empty((ref.height, ref.width), dtype="float32")

    with rasterio.open("B11.tif") as src:      # 20 m source
        reproject(
            source=rasterio.band(src, 1),
            destination=dst,
            src_transform=src.transform, src_crs=src.crs,
            dst_transform=ref.transform, dst_crs=ref.crs,
            resampling=Resampling.bilinear,
        )
# dst is now a 10 m array aligned to B04

2. Nearest-neighbour for the SCL mask

The Scene Classification Layer is a 20 m categorical band. Resample it to 10 m with Resampling.nearest so class codes stay valid — bilinear would produce fractional, meaningless codes:

import rioxarray
from rasterio.enums import Resampling

ref_10m = rioxarray.open_rasterio("B04.tif")
scl_20m = rioxarray.open_rasterio("SCL.tif")

scl_10m = scl_20m.rio.reproject_match(
    ref_10m, resampling=Resampling.nearest,   # preserve integer class codes
)

3. Aligning many bands at once from a STAC item

When bands arrive as asset hrefs from a catalog query, resolve the 10 m reference once and loop. Pair this with Optimizing Rasterio Window Reads for Memory Efficiency to avoid loading whole 10980×10980 tiles when you only need a subset:

import rioxarray
from rasterio.enums import Resampling

def load_aligned(reference_href: str, hrefs: list[str]):
    ref = rioxarray.open_rasterio(reference_href, masked=True)
    return [
        rioxarray.open_rasterio(h, masked=True)
        .rio.reproject_match(ref, resampling=Resampling.bilinear)
        for h in hrefs
    ]

Common Errors

ValueError: operands could not be broadcast together

You tried arithmetic on a 10 m array (e.g. 10980×10980) and a 20 m array (5490×5490) without aligning them. Run reproject_match on the coarse band against the 10 m reference first, then compute the index.

Index values look right but are spatially shifted

The grids were forced to the same shape without snapping origins, leaving a half-pixel offset. Use reproject_match with the native 10 m band as the template rather than resampling to a bare resolution=10, and assert transform equality afterwards.

Scene-edge pixels bloom into invalid values after bilinear

The nodata sentinel (often 0) was interpolated as if it were real reflectance. Open the coarse band with masked=True so nodata becomes NaN, or set the nodata value before resampling so GDAL excludes it from the interpolation.

Related

• Handling Pixel Resolution and Scaling — the parent section on resolution, transforms, and windowed reads across raster grids.
• Choosing the Right Resampling Method for Sentinel-2 — how to pick bilinear, cubic, or nearest for each band type.
• Advanced Resampling & Upscaling Techniques — the cross-section cluster on upscaling, downscaling, and kernel trade-offs.
• Optimizing Rasterio Window Reads for Memory Efficiency — read only the spatial subset you need before resampling large tiles.