How to Resample Satellite Bands to a Common Grid
Problem statement
Sentinel-2 delivers one scene as three grids: 10 m, 20 m and 60 m. Any analysis mixing a 10 m band with a 20 m band has to put them on one grid first, and that is a decision with a cost, not a formality.
Two costs, both measurable. Resampling method: upsampling the shortwave-infrared band from 20 m to 10 m, nearest neighbour and bilinear differ by a mean of 76.2 DN and a maximum of 8,278 DN, which changes the resulting water index by a mean of 0.017.
Grid choice: upsampling gives you 10 m spacing carrying 20 m information, so four adjacent pixels are one measurement wearing four hats.
Quick answer
Resample on read, against an explicit reference grid, with the method chosen per band type:
import rasterio
from rasterio.enums import Resampling
CLASS_BANDS = {"scl", "qa_pixel", "cloud"}
with rasterio.open(reference_path) as ref:
shape, transform, crs = ref.shape, ref.transform, ref.crs
bands = {}
for name, path in paths.items():
with rasterio.open(path) as ds:
assert ds.crs == crs, f"{name} is in {ds.crs}, reference is {crs}"
bands[name] = ds.read(
1, out_shape=shape,
resampling=Resampling.nearest if name in CLASS_BANDS
else Resampling.bilinear,
)out_shape makes GDAL resample during the read, using the dataset's own georeferencing. That is aligned by construction — unlike array-level tricks such as np.kron, which cannot see the transform at all.
Step-by-step solution
1. Choose the reference grid before you write any code
Three defensible choices, and they suit different work:
- The finest band (10 m). Keeps the visible and near-infrared detail; upsamples SWIR. Best for mapping, visual products and anything where boundaries matter.
- The coarsest band you are actually using (20 m). Downsamples the fine bands, invents nothing. Best for statistics, classification and any confidence interval computed from pixel counts.
- An external project grid. If the output must stack with other layers, match them and save a resampling step later.
The one bad choice is "whatever the first band I happened to open was".
2. Do not compute the target shape by arithmetic
The obvious approach is to double the 20 m array. It fails, because the two windows are cut against different grids:
red (10 m): (1017, 1087) origin (422200.0, 5886250.0)
swir16 (20 m): (509, 543) origin (422200.0, 5886260.0)
543 × 2 is 1086, not 1087, and the origins differ by 10 m. Take the shape and transform from the reference dataset, never from a calculation.
3. Match the method to what the band contains
| band type | method | why |
|---|---|---|
| reflectance, upsampling | bilinear |
smooth, no blockiness |
| reflectance, downsampling by ≥2 | average |
uses every source pixel |
| classification, quality, any code | nearest |
never invents a class |
| already-computed index | bilinear, reluctantly |
see below |
average matters when going coarse. bilinear downsampling by a factor of three samples the source grid and ignores most of it; average uses every contributing pixel, which is what "a 30 m pixel over this 10 m data" should mean.
4. Never interpolate class codes
Halfway between class 4 (vegetation) and class 6 (water) is class 5 (bare soil). Bilinear resampling of a classification band produces fractional values that get cast back to integers, silently relabelling every boundary pixel.
5. Compute indices at native resolution where you can
Resampling an index is not the same as computing the index from resampled bands, because the normalised difference is non-linear. The difference is small in smooth areas and largest exactly at edges — coastlines, field boundaries, cloud edges — which is where people look.
Where the bands share a resolution, compute the index first and resample the result only if you must. Where they do not, you have no choice but to resample a band, and then bilinear on the reflectance is the lesser evil.
Code examples
Example 1 — align a whole scene onto one grid
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.warp import reproject
CLASS_BANDS = {"scl", "qa_pixel", "qa", "cloud"}
def align_to_reference(paths, reference, downsample_method=Resampling.average):
"""Every band on the reference grid, method chosen per band."""
with rasterio.open(paths[reference]) as ref:
shape, transform, crs = ref.shape, ref.transform, ref.crs
ref_res = abs(ref.transform.a)
out, notes = {}, []
for name, path in paths.items():
with rasterio.open(path) as ds:
src_res = abs(ds.transform.a)
if name in CLASS_BANDS:
method = Resampling.nearest
elif src_res < ref_res / 1.5:
method = downsample_method # going coarser
else:
method = Resampling.bilinear # going finer or equal
if ds.crs != crs:
dest = np.zeros(shape, dtype=ds.dtypes[0])
reproject(rasterio.band(ds, 1), dest, dst_transform=transform,
dst_crs=crs, resampling=method)
out[name] = dest
notes.append(f"{name}: reprojected {ds.crs} -> {crs}")
else:
out[name] = ds.read(1, out_shape=shape, resampling=method)
notes.append(f"{name}: {src_res:.0f} m -> {ref_res:.0f} m "
f"({method.name})")
for note in notes:
print(" " + note)
return out, {"shape": shape, "transform": list(transform)[:6], "crs": str(crs)} red: 10 m -> 10 m (bilinear)
nir: 10 m -> 10 m (bilinear)
swir16: 20 m -> 10 m (bilinear)
scl: 20 m -> 10 m (nearest)
The src_res < ref_res / 1.5 test is what picks average for genuine downsampling and leaves everything else on bilinear. The 1.5 rather than 1.0 avoids flipping method on a rounding difference.
Example 2 — measuring what your method choice costs
import numpy as np
import rasterio
from rasterio.enums import Resampling
def resampling_cost(path, shape, methods=(Resampling.nearest,
Resampling.bilinear,
Resampling.cubic)):
"""How much do the candidate methods disagree on this band?"""
arrays = {}
with rasterio.open(path) as ds:
for method in methods:
arrays[method.name] = ds.read(1, out_shape=shape,
resampling=method).astype("float32")
names = list(arrays)
for i, a in enumerate(names):
for b in names[i + 1:]:
diff = np.abs(arrays[a] - arrays[b])
print(f" {a:9} vs {b:9}: mean {diff.mean():7.1f} "
f"p99 {np.percentile(diff, 99):8.1f} max {diff.max():8.1f}")
return arrays nearest vs bilinear : mean 76.2 p99 490.0 max 8278.0
nearest vs cubic : mean 74.5 p99 489.0 max 9320.0
bilinear vs cubic : mean 20.4 p99 114.0 max 1540.0
Run this once on a representative scene. Two readings here. Nearest disagrees with both smooth methods by about the same amount (76.2 and 74.5), while bilinear and cubic agree with each other four times more closely (20.4) — so the real decision is blocky or smooth, not which smooth kernel. And the maximum disagreement, 8,278 DN, is most of the band's range, which is what happens at a cloud edge. If the methods agree to within your analysis tolerance, stop worrying; if they do not, the choice is a documented parameter.
Example 3 — a project grid every product is written to
import json
import rasterio
from rasterio.transform import from_origin
from rasterio.warp import calculate_default_transform
def project_grid(bounds, resolution, crs, snap=True):
"""Define once, warp everything to it, and stop having alignment bugs."""
left, bottom, right, top = bounds
if snap:
left = (left // resolution) * resolution
top = ((top // resolution) + 1) * resolution
width = int((right - left) // resolution)
height = int((top - bottom) // resolution)
transform = from_origin(left, top, resolution, resolution)
grid = {"crs": str(crs), "transform": list(transform)[:6],
"width": width, "height": height, "resolution": resolution}
print(f" {width} x {height} at {resolution} m, origin ({left}, {top})")
return grid
def to_project_grid(src_path, grid, method=None):
"""Warp any raster onto the project grid."""
import numpy as np
from rasterio.warp import reproject, Resampling
from affine import Affine
transform = Affine(*grid["transform"])
with rasterio.open(src_path) as ds:
dest = np.zeros((grid["height"], grid["width"]), dtype="float32")
reproject(rasterio.band(ds, 1), dest,
dst_transform=transform, dst_crs=grid["crs"],
resampling=method or Resampling.bilinear,
src_nodata=ds.nodata, dst_nodata=np.nan)
return destSnapping the origin to a multiple of the resolution is the detail that makes grids at different resolutions share pixel corners. Without it, a 10 m and a 30 m product over the same area can be offset by a few metres forever.
Explanation
Why out_shape is better than resampling afterwards
Three reasons, in increasing order of importance.
It never materialises the full-resolution intermediate array, which for a full Sentinel-2 band is 241 MB as float32.
It uses the dataset's own transform, so the output is georeferenced correctly by construction. Array-level upsampling — np.kron, np.repeat, slicing — operates on indices and cannot correct the 10 m origin difference measured above.
And where the file has overviews, GDAL can read a lower-resolution level directly instead of reading full resolution and discarding pixels. On a COG that turns a downsampling read into a fraction of the bytes.
Why upsampling is not free information
A 20 m measurement upsampled to 10 m produces four pixels from one observation. Every value is real, and neighbouring values are not independent.
Two consequences. A classifier trained on upsampled bands sees four correlated samples where there was one, so any accuracy estimate based on pixel counts is optimistic. And a confidence interval computed from n pixels is too narrow by roughly the square root of the upsampling factor.
For anything statistical, downsample to the coarsest band you are actually using. For anything cartographic, upsample and enjoy the crisper boundaries.
Why the index-of-resampled and resampled-index differ
The normalised difference is non-linear, so f(mean(x)) ≠ mean(f(x)). Interpolating reflectances and then dividing is not the same as dividing and then interpolating.
The gap is negligible where the field is smooth and largest across sharp boundaries, because that is where the interpolation is doing the most work. Since sharp boundaries are usually the features of interest, the difference shows up exactly where it is least welcome.
Why nodata needs explicit handling
Bilinear interpolation of a pixel adjacent to nodata mixes the fill value into a real pixel. If nodata is 0, the neighbouring real values get dragged towards zero in a one-pixel band around every hole.
reproject handles this correctly when told src_nodata and dst_nodata. ds.read(out_shape=...) does not know about masking unless you use masked=True. Set nodata explicitly, or mask before resampling and accept the mask spreading by one pixel — which is usually the safer error.
Edge cases or notes
- Take the shape and transform from the reference dataset, never from arithmetic on resolutions.
np.kronand slicing cannot correct an origin offset. They work on indices, not coordinates.out_shaperesamples only within one CRS. Userasterio.warp.reprojectacross CRSs.averagefor downsampling by 2× or more; bilinear ignores most source pixels.nearestfor every class or quality band, always.- Bilinear spreads nodata by one pixel. Pass
src_nodata/dst_nodataor mask first. - Reproject once. Warping through an intermediate grid compounds interpolation error.
- Check the output range. Cubic resampling can overshoot beyond the input's minimum and maximum.
Internal links
- Satellite bands have different shapes or do not align — the errors this prevents
- How to load Sentinel-2 bands into Python as an analysis-ready array — where this fits in a loader
- 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 index that needs one grid
- Raster and vector do not line up in Python — the same problem across data models
- The raster data model explained — transforms and what they encode
FAQ
How do I resample satellite bands to the same grid?
Open a reference band, take its shape and transform, and read every other band with out_shape=reference_shape and a resampling method chosen for the band type.
Should I upsample to 10 m or downsample to 20 m?
Upsample for mapping and visual products; downsample for statistics and classification. Upsampled pixels are not independent observations.
Which resampling method should I use?
Bilinear for continuous bands going finer, average for continuous bands going coarser by 2× or more, and nearest for every classification or quality band.
How much does the resampling method matter?
Measured on a shortwave-infrared band upsampled to 10 m, nearest and bilinear differed by 76 DN on average and 8,278 DN at most, changing NDWI by an average of 0.017.
Can I just repeat each coarse pixel twice?
No. The 10 m and 20 m windows over one area were 1087 and 543 columns — doubling gives 1086 — and array repetition cannot correct a 10 m origin offset.
Should I compute an index before or after resampling?
Before, where the bands allow it. The normalised difference is non-linear, so resampling then computing gives a different answer, and the difference is largest at boundaries.
What happens to nodata when I resample?
Bilinear interpolation mixes it into neighbouring real pixels. Pass src_nodata and dst_nodata to reproject, or mask first and accept the mask growing by a pixel.