Satellite Bands Have Different Shapes or Do Not Align

Problem statement

You read two bands of the same scene, subtract them, and NumPy refuses:

import rasterio

red = rasterio.open("clear_red.tif").read(1).astype("float32")
swir = rasterio.open("clear_swir16.tif").read(1).astype("float32")
try:
    red - swir
except ValueError as e:
    print("ValueError:", e)

ValueError: operands could not be broadcast together with shapes (1017,1087) (509,543)

That is the friendly version. The dangerous version is when the shapes happen to match and the pixels still do not: two bands cut from the same bounding box, at the same resolution, whose origins differ by half a pixel. Then the arithmetic succeeds and every result is offset from reality by a few metres.

Both come from the same root cause. Bands of one satellite product are separate rasters with separate grids, and "the same area" is not the same thing as "the same pixels".

Quick answer

Never compare shapes. Compare transforms:

import rasterio


def same_grid(*paths):
    """Do these rasters describe the identical pixels?"""
    grids = []
    for path in paths:
        with rasterio.open(path) as ds:
            grids.append((path, ds.crs, ds.shape, tuple(round(v, 6)
                                                        for v in ds.transform[:6])))
    first = grids[0]
    for path, crs, shape, transform in grids:
        ok = (crs, shape, transform) == (first[1], first[2], first[3])
        print(f"  {path.split('/')[-1]:22} {str(shape):14} "
              f"origin ({transform[2]:.1f}, {transform[5]:.1f})  "
              f"{'match' if ok else 'DIFFERS'}")
    return all((c, s, t) == (first[1], first[2], first[3]) for _, c, s, t in grids)


same_grid("clear_red.tif", "clear_swir16.tif")

  clear_red.tif          (1017, 1087)   origin (422200.0, 5886250.0)  match
  clear_swir16.tif       (509, 543)     origin (422200.0, 5886260.0)  DIFFERS

Same CRS, same bounding box request, different grid and a 10 m difference in origin.

Which mismatch, no assumed factors, resampling on read

1. Identify which mismatch you have

There are four, and they need different fixes:

2. Do not assume a factor of two

The obvious fix for a 10 m/20 m mismatch is to repeat each coarse pixel:

import numpy as np

swir10 = np.kron(swir, np.ones((2, 2)))
print("red", red.shape, " swir doubled", swir10.shape)
try:
    red[swir10 > 1000]                 # "red where SWIR is bright"
except IndexError as e:
    print("IndexError:", e)

red (1017, 1087)  swir doubled (1018, 1086)
IndexError: boolean index did not match indexed array along axis 0; size of axis is 1017 but size of corresponding boolean axis is 1018

The 10 m window over this area is 1087 columns and the 20 m window is 543. Doubling gives 1086 columns — and 1018 rows against 1017, which is the axis NumPy trips on first. The bounding box was cut independently against each grid, and each rounded to its own pixel boundaries.

Even when the arithmetic works, np.kron places the upsampled pixels according to array index, not according to the georeferencing — so a half-pixel origin difference stays uncorrected.

3. Resample on read against a reference grid

import rasterio
from rasterio.enums import Resampling

with rasterio.open("clear_red.tif") as ref:
    shape, transform, crs = ref.shape, ref.transform, ref.crs

with rasterio.open("clear_swir16.tif") as ds:
    swir = ds.read(1, out_shape=shape, resampling=Resampling.bilinear)

GDAL resamples using the dataset's own transform, so the result is aligned by construction rather than by array arithmetic. It also never materialises the full-resolution intermediate.

This is correct only when both rasters cover the same ground and share a CRS. When they do not, use rasterio.warp.reproject, which takes both transforms explicitly.

4. Choose the resampling method by what the band contains

resampling = Resampling.nearest if band == "scl" else Resampling.bilinear

Class codes must use nearest neighbour. Interpolating between class 4 (vegetation) and class 6 (water) gives class 5 (bare soil) — a class that was never observed at that location.

For continuous bands the choice costs real accuracy. Measured on the shortwave-infrared band upsampled from 20 m to 10 m:

with rasterio.open("clear_swir16.tif") as ds:
    near, bil = (ds.read(1, out_shape=shape, resampling=m).astype("float32") / 1e4
                 for m in (Resampling.nearest, Resampling.bilinear))
nir = rasterio.open("clear_nir.tif").read(1).astype("float32") / 1e4
ndwi = lambda sw: (nir - sw) / (nir + sw)
print(f"nearest vs bilinear: mean difference {np.abs(near - bil).mean() * 1e4:.1f} DN, "
      f"max {np.abs(near - bil).max() * 1e4:,.0f} DN")
print(f"resulting index difference (NDWI): mean {np.abs(ndwi(near) - ndwi(bil)).mean():.4f}")

nearest vs bilinear: mean difference 41.8 DN, max 1,853 DN
resulting index difference (NDWI): mean 0.0088

An index difference of 0.0088 is small against a class threshold and large against an inter-annual trend.

5. Verify after resampling, not before

assert swir.shape == shape
with rasterio.open("clear_swir16.tif") as ds:
    assert ds.crs == crs

The shape assertion is trivially true after out_shape. The CRS assertion is the one that catches a scene straddling two UTM zones.

Aligning bands onto one reference

The data used on this page

The least cloudy Sentinel-2 L2A scene of tile 30UVD over Snowdon, June–September 2025, from Earth Search: the setup cuts the red, near-infrared, shortwave-infrared (B11) and scene classification bands to the area 4.16–4.00° W, 53.03–53.12° N, each at its native resolution and against its own grid, and writes them as clear_<band>.tif.

import numpy as np
import rasterio
from pystac_client import Client
from rasterio.warp import transform_bounds

AOI = (-4.16, 53.03, -4.00, 53.12)
items = sorted(Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}}).items(), key=lambda i: i.properties["eo:cloud_cover"])
item = items[0]
print(item.id, f"({item.properties['eo:cloud_cover']:.1f}% cloud)")

paths = {}
for name in ("red", "nir", "swir16", "scl"):
    with rasterio.open(item.assets[name].href) as ds:
        win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                           transform=ds.transform).round_offsets().round_lengths()
        data = ds.read(1, window=win)
        profile = {"driver": "GTiff", "dtype": data.dtype, "count": 1, "crs": ds.crs, "nodata": ds.nodata,
                   "width": data.shape[1], "height": data.shape[0], "transform": ds.window_transform(win)}
    paths[name] = f"clear_{name}.tif"
    with rasterio.open(paths[name], "w", **profile) as out:
        out.write(data, 1)
    t = profile["transform"]
    print(f"{name:7} ({t.a:.0f} m): {data.shape}   origin ({t.c}, {t.f})")
reference_path = paths["red"]

S2C_30UVD_20250712_0_L2A (0.1% cloud)
red     (10 m): (1017, 1087)   origin (422200.0, 5886250.0)
nir     (10 m): (1017, 1087)   origin (422200.0, 5886250.0)
swir16  (20 m): (509, 543)   origin (422200.0, 5886260.0)
scl     (20 m): (509, 543)   origin (422200.0, 5886260.0)

Example 1 — align any set of bands onto one reference

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

CLASS_BANDS = {"scl", "qa", "cloud", "classification"}


def align_bands(paths, reference, method=None):
    """Every band on the reference band's grid, reprojecting if the CRS differs."""
    with rasterio.open(paths[reference]) as ref:
        shape, transform, crs = ref.shape, ref.transform, ref.crs

    out = {}
    for name, path in paths.items():
        with rasterio.open(path) as ds:
            resampling = method or (Resampling.nearest if name in CLASS_BANDS
                                    else Resampling.bilinear)
            if ds.crs == crs:
                # same CRS: resample on read, aligned by the dataset transform
                out[name] = ds.read(1, out_shape=shape, resampling=resampling)
            else:
                # different CRS: a full reprojection is required
                dest = np.zeros(shape, dtype=ds.dtypes[0])
                reproject(rasterio.band(ds, 1), dest,
                          dst_transform=transform, dst_crs=crs,
                          resampling=resampling)
                out[name] = dest
                print(f"  {name}: reprojected from {ds.crs} to {crs}")

        print(f"  {name:10} -> {out[name].shape} "
              f"({resampling.name})")
    return out, {"shape": shape, "transform": list(transform)[:6], "crs": str(crs)}


aligned, grid = align_bands(paths, "red")

  red        -> (1017, 1087) (bilinear)
  nir        -> (1017, 1087) (bilinear)
  swir16     -> (1017, 1087) (bilinear)
  scl        -> (1017, 1087) (nearest)

Reading out_shape for the reference band itself is a no-op, so the reference can stay in the dictionary and the code has no special case.

Example 2 — detect a silent misalignment

import numpy as np
import rasterio


def alignment_offset(path_a, path_b):
    """How far apart are two grids, in pixels of the first?"""
    with rasterio.open(path_a) as a, rasterio.open(path_b) as b:
        if a.crs != b.crs:
            return {"error": f"different CRS: {a.crs} vs {b.crs}"}

        dx = (b.transform.c - a.transform.c) / a.transform.a
        dy = (b.transform.f - a.transform.f) / a.transform.e
        scale_x = b.transform.a / a.transform.a
        scale_y = b.transform.e / a.transform.e

    result = {
        "offset_px": (round(dx, 3), round(dy, 3)),
        "offset_m": (round(dx * a_res, 1) if (a_res := abs(a.transform.a)) else None,
                     None),
        "pixel_ratio": (round(scale_x, 3), round(scale_y, 3)),
        "aligned": abs(dx % 1) < 1e-6 and abs(dy % 1) < 1e-6,
    }
    print(f"  offset {result['offset_px']} px, "
          f"pixel ratio {result['pixel_ratio']}, "
          f"aligned on the grid: {result['aligned']}")
    return result


alignment_offset("clear_red.tif", "clear_swir16.tif")

  offset (0.0, -1.0) px, pixel ratio (2.0, 2.0), aligned on the grid: True

A whole-pixel offset is recoverable. A fractional one means one of the two rasters was written from a warp that did not snap to the grid, and it will smear every subsequent resampling.

Example 3 — a reusable reference grid for a whole project

import json
import rasterio
from rasterio.transform import from_origin


def define_grid(bounds, resolution, crs, snap=True):
    """One grid definition every product in the project is written to."""
    left, bottom, right, top = bounds
    if snap:
        # snap the origin to a multiple of the resolution so any two
        # products at any resolution share pixel corners
        left = (left // resolution) * resolution
        top = ((top // resolution) + 1) * resolution

    width = int((right - left) // resolution)
    height = int((top - bottom) // resolution)

    grid = {
        "crs": str(crs),
        "transform": list(from_origin(left, top, resolution, resolution))[:6],
        "width": width, "height": height, "resolution": resolution,
    }
    print(f"  {width} x {height} at {resolution} m, origin ({left}, {top})")
    return grid


with rasterio.open("clear_red.tif") as ds:
    project = define_grid(ds.bounds, 10, ds.crs)

  1087 x 1018 at 10 m, origin (422200.0, 5886260.0)

Defining the grid once, at the top of the project, and warping everything to it removes this entire class of problem. It costs one resampling of each input and buys arithmetic that is correct by construction rather than by inspection.

Example 4 — where the off-by-one comes from

from rasterio.warp import transform_bounds

with rasterio.open(item.assets["red"].href) as ten, rasterio.open(item.assets["swir16"].href) as twenty:
    left, bottom, right, top = transform_bounds("EPSG:4326", ten.crs, *AOI)
    print(f"requested box: {right - left:,.1f} m wide, {top - bottom:,.1f} m tall")
    for name, ds in (("10 m", ten), ("20 m", twenty)):
        w = rasterio.windows.from_bounds(left, bottom, right, top, transform=ds.transform)
        r = w.round_offsets().round_lengths()
        print(f"  {name}: {w.height:.2f} x {w.width:.2f} pixels -> {r.height} x {r.width}, "
              f"top edge {top:,.1f} snapped to {ds.window_transform(r).f:,.1f}")

native = 4 * 1 + 6 * (10 / 20) ** 2 + 3 * (10 / 60) ** 2      # band-grids at 10 m equivalent
print(f"13 bands: {native:.2f} 10 m grids' worth stored natively, 13 if all were resampled to 10 m "
      f"({13 / native:.1f}x)")

requested box: 10,869.4 m wide, 10,172.9 m tall
  10 m: 1017.29 x 1086.94 pixels -> 1017 x 1087, top edge 5,886,248.3 snapped to 5,886,250.0
  20 m: 508.65 x 543.47 pixels -> 509 x 543, top edge 5,886,248.3 snapped to 5,886,260.0
13 bands: 5.58 10 m grids' worth stored natively, 13 if all were resampled to 10 m (2.3x)

Why bands are separate grids

Why bands are separate grids in the first place

Sentinel-2's bands are acquired by different detector arrays with different pixel sizes: 10 m for four bands, 20 m for six, 60 m for three. Storing them on one grid would mean either upsampling the coarse bands — inventing detail and storing 2.3 times as much (Example 4) — or downsampling the fine ones and throwing away what makes them useful.

Keeping them separate defers the decision to whoever knows the analysis. That is the right call, and it means every user makes the decision, including the users who did not realise there was one.

Why the shapes come out mismatched by one

rasterio.windows.from_bounds computes a window in fractional pixels, and .round_lengths() rounds to whole ones. A 10,173 m span is 1017.3 pixels at 10 m and 508.65 at 20 m, rounding to 1017 and 509. Doubling 509 gives 1018, not 1017.

There is no rounding rule that avoids this in general, because the two grids genuinely do not have a common pixel boundary at every requested bound. The fix is not better rounding; it is to define one grid and resample onto it.

Why the origin differs by 10 m

Each window snaps to its own grid's pixel edges. The 10 m grid has an edge every 10 m; the 20 m grid every 20 m. .round_offsets() floors the pixel offset, so each window's top edge moves up to the grid line at or above the request. Here the box's top at 5,886,248.3 m became 5,886,250 on the 10 m grid and 5,886,260 on the 20 m grid, 11.7 m above what was asked for (Example 4).

This is invisible in the shape and visible in the transform, which is why the check is on the transform.

Why nearest neighbour on class bands is not a preference

Resampling interpolates. For reflectance that is meaningful — halfway between 0.2 and 0.4 is 0.3, a plausible reflectance. For a classification band the codes are labels, and their numeric order carries no meaning: 4 is vegetation, 5 is bare soil, 6 is water, and the average of vegetation and water is not bare soil.

Bilinear resampling of a class band produces non-integer values that then get cast back to integers, silently relabelling boundary pixels as whatever class number happens to lie between them.

Shapes that match but do not align

• Matching shapes do not prove alignment. Compare transform and crs.
• np.kron and array slicing ignore georeferencing. They cannot correct an origin offset.
• out_shape only works within one CRS. Across CRSs use rasterio.warp.reproject.
• Nearest for classes, bilinear or cubic for continuous, average for downsampling by a large factor.
• Round windows with .round_offsets().round_lengths(), and expect off-by-one across resolutions anyway.
• A scene can straddle two UTM zones; the same date can arrive as two items in different CRSs.
• Reprojecting twice compounds error. Warp each input once, from its native grid to the project grid.
• Check nodata survives resampling. Bilinear interpolation of a nodata edge spreads the fill value into real pixels.

Internal links

• How to resample satellite bands to a common grid — the full procedure and its costs
• How to load Sentinel-2 bands into Python as an analysis-ready array — a loader that aligns by construction
• Raster and vector do not line up in Python — the same problem across data models
• Raster resampling explained — what each method does to the values
• How to reproject a raster in Python with rasterio — when the CRSs differ
• Spectral bands explained — why the resolutions differ at all
• How to calculate NDVI from Sentinel-2 in Python — the first calculation this breaks
• The raster data model explained — what a transform is

FAQ

Why do Sentinel-2 bands have different shapes?

Because they have different native resolutions — 10 m, 20 m and 60 m — and each band is stored on its own grid rather than resampled to a common one.

Can I just repeat each 20 m pixel twice?

No. The two windows are cut against different grids, so the counts often do not divide exactly — 543 × 2 is 1086 against a 10 m width of 1087 — and array-index repetition cannot correct an origin offset.

Why do my bands have the same shape but give shifted results?

Different origins. Two windows cut from the same bounding box snap to different pixel edges. Compare ds.transform, not ds.shape.

Which resampling method should I use?

Nearest for classification and quality bands, bilinear for continuous bands when upsampling, average when downsampling by a large factor.

How much does the resampling method matter?

Measured on the shortwave-infrared band upsampled to 10 m, nearest and bilinear differed by 42 DN on average and produced NDWI values differing by 0.0088 on average.

Should I upsample the coarse bands or downsample the fine ones?

Upsample for mapping, downsample for statistics. Upsampling gives fine spacing with coarse information, which makes neighbouring pixels correlated.

What if the two bands are in different CRSs?

Use rasterio.warp.reproject with both transforms and CRSs given explicitly. out_shape on read only resamples within one CRS.