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 41.8 DN and a maximum of 1,853 DN on an almost cloud-free Snowdon scene, which changes the resulting water index by a mean of 0.0088 and by 0.068 at the 99th percentile.

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.

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)

Choose the reference grid, no shape arithmetic

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

average matters when going coarse, though not for the reason usually given. GDAL widens the bilinear kernel when it downsamples, so bilinear does read every source pixel: one bright pixel placed in each of the nine positions of a 3 × 3 block reached the 30 m output all nine times (Example 4). But bilinear also spreads each pixel over 2.8 output cells on average, and cubic over 8.6, while average keeps it in the one cell it falls in — which is what "a 30 m pixel over this 10 m data" should mean. On the red band at 30 m, bilinear differed from average by 11.3 DN on average; nearest, which reached only one of the nine positions, by 33.9.

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. Taking 10 m NDVI to 30 m, index-first and bands-first differed by 0.0003 in the smoothest half of the cells and by 0.0183 in the most varied tenth (Example 4).

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.

Aligning a whole scene and measuring the cost

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 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)}


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

  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


arrays = resampling_cost(paths["swir16"], grid["shape"])
nir = aligned["nir"] / 1e4


def ndwi(swir):
    """Gao's water index from NIR and SWIR reflectance."""
    return (nir - swir / 1e4) / (nir + swir / 1e4)


index_diff = np.abs(ndwi(arrays["nearest"]) - ndwi(arrays["bilinear"]))
print(f"  NDWI, nearest vs bilinear SWIR: mean {index_diff.mean():.4f}  p99 {np.percentile(index_diff, 99):.4f}")

  nearest   vs bilinear : mean    41.8  p99    275.0  max   1853.0
  nearest   vs cubic    : mean    40.2  p99    271.0  max   1916.0
  bilinear  vs cubic    : mean    13.1  p99     74.0  max    390.0
  NDWI, nearest vs bilinear SWIR: mean 0.0088  p99 0.0684

Run this once on a representative scene. Two readings here. Nearest disagrees with both smooth methods by about the same amount (41.8 and 40.2), while bilinear and cubic agree with each other three times more closely (13.1) — so the real decision is blocky or smooth, not which smooth kernel. And the maximum disagreement, 1,853 DN, is about a quarter of the band's 77–7,016 DN range in this window; it sits where a dark area meets bare soil on the mountain (Example 4). This scene is 0.1% cloud, and the SCL band finds none in the window. 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 dest


with rasterio.open(paths["red"]) as ds:
    red_bounds, red_crs = ds.bounds, ds.crs
grid10 = project_grid(red_bounds, 10, red_crs)
grid30 = project_grid(red_bounds, 30, red_crs)
swir30 = to_project_grid(paths["swir16"], grid30, Resampling.average)
print(f"  swir16 on the 30 m grid: {swir30.shape}, {np.isnan(swir30).sum()} empty cells")

  1087 x 1018 at 10 m, origin (422200.0, 5886260.0)
  362 x 339 at 30 m, origin (422190.0, 5886270.0)
  swir16 on the 30 m grid: (339, 362), 0 empty cells

Snapping 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.

Example 4 — what GDAL's methods actually do

from rasterio.io import MemoryFile
from rasterio.transform import from_origin


def in_memory(array, transform, nodata=None):
    profile = {"driver": "GTiff", "dtype": array.dtype, "count": 1, "width": array.shape[1],
               "height": array.shape[0], "crs": "EPSG:32630", "transform": transform, "nodata": nodata}
    memfile = MemoryFile()
    with memfile.open(**profile) as dst:
        dst.write(array, 1)
    return memfile


# 1. Downsampling 10 m -> 30 m: one bright pixel in each of the nine positions of a 3 x 3 block
impulses = np.zeros((90, 270), "float32")
for k in range(9):
    impulses[30 + k // 3, 30 * k + 3 + k % 3] = 1000
with in_memory(impulses, from_origin(0, 900, 10, 10)).open() as ds:
    for method in (Resampling.nearest, Resampling.bilinear, Resampling.cubic, Resampling.average):
        out = ds.read(1, out_shape=(30, 90), resampling=method)
        seen = sum(out[9:12, 10 * k:10 * k + 3].any() for k in range(9))
        print(f"downsample x3, {method.name:8}: {seen}/9 source positions reach the output, "
              f"each spreads over {np.count_nonzero(out) / max(seen, 1):.1f} cells")

with rasterio.open(paths["red"]) as ds:
    thirty = {m.name: ds.read(1, out_shape=(ds.height // 3, ds.width // 3), resampling=m).astype("float32")
              for m in (Resampling.nearest, Resampling.bilinear, Resampling.average)}
for m in ("nearest", "bilinear"):
    print(f"red at 30 m, {m} vs average: mean difference {np.abs(thirty[m] - thirty['average']).mean():.1f} DN")

# 2. Nodata: a 1,000 DN field with a hole of zeros, upsampled 20 m -> 10 m
hole = np.full((20, 20), 1000, "uint16")
hole[8:12, 8:12] = 0
for nodata in (0, None):
    with in_memory(hole, from_origin(0, 400, 20, 20), nodata).open() as ds:
        up = ds.read(1, out_shape=(40, 40), resampling=Resampling.bilinear)
    real = up[up > 0]
    print(f"nodata={nodata}: lowest real pixel {real.min()}, {np.count_nonzero(real < 1000)} dragged below 1,000")
for kwargs in ({}, {"src_nodata": 0, "dst_nodata": 0}):
    up = np.zeros((40, 40), "uint16")
    reproject(hole, up, src_transform=from_origin(0, 400, 20, 20), src_crs="EPSG:32630",
              dst_transform=from_origin(0, 400, 10, 10), dst_crs="EPSG:32630", resampling=Resampling.bilinear, **kwargs)
    real = up[up > 0]
    print(f"reproject {kwargs or 'without src_nodata'}: lowest real pixel {real.min()}, "
          f"{np.count_nonzero(real < 1000)} dragged below 1,000")

# 3. Cubic overshoot on the SWIR band
with rasterio.open(paths["swir16"]) as ds:
    native = ds.read(1)
print(f"swir16 native {native.min()}-{native.max()} DN, cubic to 10 m "
      f"{arrays['cubic'].min():.0f}-{arrays['cubic'].max():.0f}, bilinear "
      f"{arrays['bilinear'].min():.0f}-{arrays['bilinear'].max():.0f}")

# 4. Where nearest and bilinear disagree most
iy, ix = np.unravel_index(np.abs(arrays["nearest"] - arrays["bilinear"]).argmax(), grid["shape"])
scl_names = {2: "dark area", 3: "cloud shadow", 4: "vegetation", 5: "bare soil", 6: "water", 7: "unclassified",
             8: "cloud medium", 9: "cloud high", 10: "thin cirrus"}
around = np.unique(aligned["scl"][iy - 2:iy + 3, ix - 2:ix + 3])
print(f"largest difference at row {iy}, col {ix}: SCL classes around it: {[scl_names[c] for c in around]}")
print("SCL classes in the window:", {scl_names[c]: int(n) for c, n in zip(*np.unique(aligned["scl"], return_counts=True))})

# 5. Index then resample, against resample then index, 10 m -> 30 m
red10, nir10 = aligned["red"][:1017 // 3 * 3, :1087 // 3 * 3] / 1e4, aligned["nir"][:1017 // 3 * 3, :1087 // 3 * 3] / 1e4
blocks = lambda a: a.reshape(a.shape[0] // 3, 3, a.shape[1] // 3, 3)
ndvi10 = (nir10 - red10) / (nir10 + red10)
index_first = blocks(ndvi10).mean(axis=(1, 3))
bands_first = (blocks(nir10).mean(axis=(1, 3)) - blocks(red10).mean(axis=(1, 3))) / \
              (blocks(nir10).mean(axis=(1, 3)) + blocks(red10).mean(axis=(1, 3)))
gap, spread = np.abs(index_first - bands_first), blocks(ndvi10).std(axis=(1, 3))
edge, flat = spread >= np.percentile(spread, 90), spread <= np.percentile(spread, 50)
print(f"NDVI index-first vs bands-first at 30 m: mean gap {gap.mean():.4f}; "
      f"{gap[edge].mean():.4f} in the 10% most varied cells, {gap[flat].mean():.4f} in the smoothest half")

# 6. A full Sentinel-2 tile at 10 m
print(f"one band of a 10,980 x 10,980 tile as float32: {10980 * 10980 * 4 / 1e6:.0f} MB")

downsample x3, nearest : 1/9 source positions reach the output, each spreads over 1.0 cells
downsample x3, bilinear: 9/9 source positions reach the output, each spreads over 2.8 cells
downsample x3, cubic   : 9/9 source positions reach the output, each spreads over 8.6 cells
downsample x3, average : 9/9 source positions reach the output, each spreads over 1.0 cells
red at 30 m, nearest vs average: mean difference 33.9 DN
red at 30 m, bilinear vs average: mean difference 11.3 DN
nodata=0: lowest real pixel 1000, 0 dragged below 1,000
nodata=None: lowest real pixel 250, 64 dragged below 1,000
reproject without src_nodata: lowest real pixel 250, 64 dragged below 1,000
reproject {'src_nodata': 0, 'dst_nodata': 0}: lowest real pixel 1000, 0 dragged below 1,000
swir16 native 77-7016 DN, cubic to 10 m 48-7008, bilinear 79-6962
largest difference at row 475, col 630: SCL classes around it: ['dark area', 'bare soil']
SCL classes in the window: {'dark area': 6868, 'vegetation': 1029513, 'bare soil': 42308, 'water': 23322, 'unclassified': 3468}
NDVI index-first vs bands-first at 30 m: mean gap 0.0027; 0.0183 in the 10% most varied cells, 0.0003 in the smoothest half
one band of a 10,980 x 10,980 tile as float32: 482 MB

Why out_shape beats resampling afterwards

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 at 10 m is 482 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 a fill value mixes the fill into a real pixel — if GDAL does not know it is a fill. In Example 4, a 1,000 DN field with a 4 × 4 hole of zeros, upsampled 20 m to 10 m, kept every real pixel at 1,000 when the file declared nodata=0; with no nodata declared, 64 pixels around the hole were dragged down, the lowest to 250.

So ds.read(out_shape=...) honours the nodata value the file declares, and reproject honours the one you pass as src_nodata (and writes dst_nodata). The failure is a file that uses 0 or −9999 as fill without saying so. Check ds.nodata before resampling; if it is None and the file has holes, set it or mask first.

Reference transforms and origin offsets

• Take the shape and transform from the reference dataset, never from arithmetic on resolutions.
• np.kron and slicing cannot correct an origin offset. They work on indices, not coordinates.
• out_shape resamples only within one CRS. Use rasterio.warp.reproject across CRSs.
• average for downsampling by 2× or more; GDAL's bilinear reads every source pixel but smears each across neighbouring output cells.
• nearest for every class or quality band, always.
• Bilinear spreads an undeclared fill value. Check ds.nodata; pass src_nodata/dst_nodata to reproject.
• 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 42 DN on average and 1,853 DN at most, changing NDWI by an average of 0.0088.

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?

If the file declares its nodata value, GDAL leaves it out of the interpolation; if it does not, bilinear mixes the fill into neighbouring real pixels — 64 of them around a 4 × 4 hole in Example 4. Check ds.nodata, and pass src_nodata and dst_nodata to reproject.