Satellite Bands Have Different Shapes or Do Not Align
Problem statement
You read two bands of the same scene, subtract them, and NumPy refuses:
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) 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.
Step-by-step solution
1. Identify which mismatch you have
There are four, and they need different fixes:
| symptom | cause | fix |
|---|---|---|
| shapes differ by ~2× | mixed native resolutions | resample on read with out_shape |
| shapes differ by 1 | window rounding | reproject to a shared reference grid |
| shapes match, results shifted | different origins | check transform, not shape |
ValueError about CRS |
different UTM zones | reproject one band |
2. Do not assume a factor of two
The obvious fix for a 10 m/20 m mismatch is to repeat each coarse pixel:
swir10 = np.kron(swir, np.ones((2, 2)))IndexError: boolean index did not match indexed array along axis 1;
size of axis is 1087 but size of corresponding boolean axis is 1086
The 10 m window over this area is 1087 columns and the 20 m window is 543. Doubling gives 1086. 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.bilinearClass 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:
nearest vs bilinear: mean difference 76.2 DN, max 8,278 DN
resulting index difference (NDWI): mean 0.017
An index difference of 0.017 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 == crsThe shape assertion is trivially true after out_shape. The CRS assertion is the one that catches a scene straddling two UTM zones.
Code examples
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)} 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 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 1088 x 1018 at 10 m, origin (422200, 5886260)
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.
Explanation
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 tripling the storage — 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. A requested top edge that falls between two 20 m lines snaps to the nearer one, up to 10 m away from where the 10 m grid snapped.
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.
Edge cases or notes
- Matching shapes do not prove alignment. Compare
transformandcrs. np.kronand array slicing ignore georeferencing. They cannot correct an origin offset.out_shapeonly works within one CRS. Across CRSs userasterio.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
nodatasurvives 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 76 DN on average and produced NDWI values differing by 0.017 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.