Raster Resampling Explained: Nearest, Bilinear and When Each Is Wrong

Problem statement

You reproject a land-cover raster into British National Grid. ESA WorldCover over the South Pennines — Halifax, Huddersfield and the moors between them — codes each 10 m cell as one of a handful of classes: 10 tree cover, 30 grassland, 50 built-up, 80 water and so on. Reproject it twice, once with each of the two methods people reach for:

from rasterio.warp import calculate_default_transform, reproject, Resampling

def to_bng(path, method, resolution=10):
    with rasterio.open(path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, "EPSG:27700", src.width, src.height, *src.bounds, resolution=resolution)
        out = np.zeros((height, width), dtype=src.dtypes[0])
        reproject(rasterio.band(src, 1), out, dst_transform=transform, dst_crs="EPSG:27700",
                  src_nodata=src.nodata, dst_nodata=src.nodata, resampling=method)
    return out

for method in (Resampling.nearest, Resampling.bilinear):
    out = to_bng("landcover.tif", method)
    valid = out != 0
    invented = valid & ~np.isin(out, list(CLASSES))
    print(f"{method.name:8} {len(np.unique(out[valid])):>3} distinct values, "
          f"{invented.sum() / valid.sum():.1%} of cells hold a code that is not a class")
print(f"cells coded 45, halfway between cropland (40) and built-up (50): {(out == 45).sum():,}")

nearest    8 distinct values, 0.0% of cells hold a code that is not a class
bilinear  81 distinct values, 18.1% of cells hold a code that is not a class
cells coded 45, halfway between cropland (40) and built-up (50): 28,896

Nothing errored. The nearest-neighbour copy still holds the eight classes the source has. The bilinear one holds dozens of codes that WorldCover does not define: it averaged category codes, and a cell coded 45 is the arithmetic mean of a field and a street — not a class at all.

Every operation that changes a raster's grid — reprojection, clipping to a different alignment, changing resolution, overviews, mosaicking with mismatched grids — must decide what value each new cell gets. That decision is resampling. Rasterio's reproject defaults to nearest, so this is rarely an accident of defaults; it is bilinear chosen because it looks smoother, or copied from a snippet written for elevation.

Quick answer

from rasterio.enums import Resampling

The last row matters. Nearest neighbour never invents a value that was not in the source. It is the only method that is safe when you do not yet know what you are looking at.

from rasterio.warp import reproject, Resampling

reproject(
    source=src.read(1), destination=dest,
    src_transform=src.transform, src_crs=src.crs, src_nodata=src.nodata,
    dst_transform=dst_transform, dst_crs="EPSG:27700", dst_nodata=src.nodata,
    resampling=Resampling.nearest,     # ← never leave this to chance
)

What your numbers mean, what each method does, the match

1. Work out what your numbers mean

This is the whole decision, and it is about the data, not the algorithm.

for path in ("landcover.tif", "elevation.tif"):
    with rasterio.open(path) as src:
        data = src.read(1, masked=True)
    vals = np.unique(data.compressed())
    print(f"{path:14} {src.dtypes[0]:8} {len(vals):>8,} distinct values   first few: {vals[:5]}")

landcover.tif  uint8           8 distinct values   first few: [10 20 30 40 50]
elevation.tif  float32   841,820 distinct values   first few: [29.       29.5      30.       30.163267 30.249184]

Eight distinct integers in a uint8 raster is categorical, near-certainly — and WorldCover's documentation confirms each code is a label. Hundreds of thousands of distinct float32 heights is continuous. The rule of thumb is reliable: few distinct values in an integer dtype means the numbers are labels, and arithmetic on labels is meaningless.

The awkward middle case is an integer continuous raster — elevation in whole metres has hundreds or thousands of distinct int16 values. It is continuous. Count first, then decide.

2. Know what each method actually does

Nearest neighbour copies the value of the single closest source pixel.

• Output values are always a subset of input values.
• Preserves categories exactly.
• Produces visible blockiness and stair-stepped edges when upsampling.
• Loses small features when downsampling — a one-pixel stream disappears wherever no output cell centre lands on it.

Bilinear takes a distance-weighted average of the four surrounding pixels.

• Smooth, continuous output.
• Invents values between the inputs — correct for elevation, nonsense for class codes.
• Slightly blurs sharp real edges, such as a cliff or a coastline.

Cubic and Lanczos fit a curve through a 4 × 4 and a 6 × 6 neighbourhood; cubic spline fits a smoothing B-spline through 4 × 4. The difference matters on bounded data. Here is the tree cover of the same area as a percentage per 100 m cell — a quantity that cannot leave 0–100 — upsampled back to 10 m with each:

import time
from affine import Affine

cover = to_bng("landcover.tif", Resampling.nearest)
h, w = (s // 10 * 10 for s in cover.shape)
tree_pct = (cover[:h, :w] == 10).reshape(h // 10, 10, w // 10, 10).mean(axis=(1, 3)) * 100
print(f"tree cover per 100 m cell: {tree_pct.min():.0f} … {tree_pct.max():.0f}%")

for method in (Resampling.bilinear, Resampling.cubic_spline, Resampling.cubic, Resampling.lanczos):
    fine = np.empty((h, w), dtype="float32")
    start = time.perf_counter()
    reproject(tree_pct.astype("float32"), fine,
              src_transform=Affine(100, 0, 0, 0, -100, 0), src_crs="EPSG:27700",
              dst_transform=Affine(10, 0, 0, 0, -10, 0), dst_crs="EPSG:27700", resampling=method)
    outside = ((fine < 0) | (fine > 100)).mean()
    print(f"{method.name:13} {fine.min():7.1f} … {fine.max():6.1f}%   "
          f"{outside:6.1%} of cells outside 0–100   {time.perf_counter() - start:.2f} s")

tree cover per 100 m cell: 0 … 100%
bilinear          0.0 …  100.0%     0.0% of cells outside 0–100   0.12 s
cubic_spline      0.0 …  100.0%     0.0% of cells outside 0–100   0.42 s
cubic           -14.6 …  117.4%    14.1% of cells outside 0–100   0.12 s
lanczos         -28.4 …  129.6%    16.3% of cells outside 0–100   0.63 s

Bilinear and cubic spline stay inside the range of their inputs. Cubic and Lanczos overshoot: near a sharp step — a wood edge, where cover goes from 0 to 100 in one cell — they produce negative tree cover and cover above 100%, on a large share of cells. On elevation that ringing is harmless. On a percentage, a probability or an index bounded to −1…1 it is physically impossible, and code that assumed the bound will break or, worse, will not. The cost difference between these methods is small next to that.

Average, mode and sum aggregate, and are only meaningful when several source pixels fall inside one output pixel:

Resampling.average   # mean of the contributing pixels — continuous data
Resampling.mode      # most common contributing value — a map of categorical data
Resampling.sum       # total of the contributing pixels — counts

When upsampling there is one source pixel per output pixel, nothing to aggregate, and they return exactly what nearest does — checked on the tree-cover raster above, the arrays are identical.

3. Match the method to the direction

The direction of the scale change matters as much as the data type.

Upsampling (finer output) means one input pixel spread across many output pixels. No method adds information; the choice is only about how the interpolation looks. nearest gives blocks, bilinear gives smooth ramps.

Downsampling (coarser output) means many input pixels per output pixel. Taking the 10 m land cover to 100 m, each output cell has a hundred candidates, and the three reasonable methods answer different questions:

from rasterio.io import MemoryFile


def coarsen(array, factor, method, nodata=0):
    with MemoryFile() as mem:
        with mem.open(driver="GTiff", width=array.shape[1], height=array.shape[0], count=1,
                      dtype=array.dtype, crs="EPSG:27700", transform=Affine(10, 0, 0, 0, -10, 0),

                      nodata=nodata) as ds:
            ds.write(array, 1)
            return ds.read(1, out_shape=(array.shape[0] // factor, array.shape[1] // factor),
                           resampling=method)

def shares(array):
    valid = array[array != 0]
    return {code: (valid == code).mean() for code in CLASSES}

source = shares(cover)
near, mode = shares(coarsen(cover, 10, Resampling.nearest)), shares(coarsen(cover, 10, Resampling.mode))
print(f"{'class':12} {'10 m':>7} {'nearest':>8} {'mode':>7}")
for code in (10, 30, 40, 50, 60, 80, 20, 90):
    print(f"{CLASSES[code]:12} {source[code]:7.2%} {near[code]:8.2%} {mode[code]:7.2%}")


is_tree = np.where(cover == 0, np.nan, cover == 10).astype("float32")   # 1 tree, 0 other, NaN outside
tree_fraction = coarsen(is_tree, 10, Resampling.average, nodata=np.nan)
print(f"tree cover from 100 m fractions: {np.nanmean(tree_fraction):.2%}")

class           10 m  nearest    mode
tree cover    27.95%   28.01%  25.40%
grassland     46.68%   46.58%  48.20%
cropland       2.73%    2.75%   2.74%
built-up      22.03%   22.02%  23.16%
bare           0.09%    0.11%   0.06%
water          0.50%    0.50%   0.45%
shrubland      0.01%    0.01%   0.00%
wetland        0.01%    0.01%   0.00%
tree cover from 100 m fractions: 27.97%

nearest keeps one pixel in a hundred, but it keeps an unbiased sample: every class's share comes through almost unchanged, rare ones included. It is noisy cell by cell — which pixel wins depends on where the output grid falls — but not biased overall.

mode is biased, and in a predictable direction. Each output cell takes its majority class, so the dominant classes grow, the minority ones shrink, and a class that is never the local majority — shrubland, wetland — vanishes completely. That is the right behaviour for a map of what each 100 m cell mostly is, and the wrong one for any area statistic. For areas, keep the fractions: average of a 0/1 layer per class gives each 100 m cell its tree fraction, and those fractions reproduce the 10 m share to within 0.02 percentage points.

For continuous data, average uses all hundred pixels and is the right answer; sum is the right answer for counts such as population, where an average would divide the total by a hundred.

4. Set it explicitly, every time

def reproject_raster(src_path, dst_path, dst_crs, resampling, resolution, declare_nodata=True):
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds, resolution=resolution)
        nodata = src.nodata
        if nodata is None and declare_nodata and np.issubdtype(np.dtype(src.dtypes[0]), np.floating):
            nodata = np.nan       # the Copernicus DEM declares none
        profile = src.profile | {"crs": dst_crs, "transform": transform, "width": width,
                                 "height": height, "nodata": nodata}
        with rasterio.open(dst_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                reproject(rasterio.band(src, i), rasterio.band(dst, i),
                          src_nodata=src.nodata, dst_nodata=nodata, resampling=resampling)

reproject_raster("landcover.tif", "landcover_bng.tif", "EPSG:27700", Resampling.nearest, 10)
reproject_raster("landcover.tif", "landcover_bilinear.tif", "EPSG:27700", Resampling.bilinear, 10)
reproject_raster("elevation.tif", "elevation_bng.tif", "EPSG:27700", Resampling.bilinear, 30)
reproject_raster("elevation.tif", "elevation_cubic.tif", "EPSG:27700", Resampling.cubic, 30)
reproject_raster("elevation.tif", "elevation_no_nodata.tif", "EPSG:27700", Resampling.bilinear, 30,
                 declare_nodata=False)

Different layers, different methods, and the difference is not cosmetic. The nodata lines matter as much as resampling: a square of lon/lat cells does not stay square in British National Grid, and the corners of the new grid that fall outside the source have to hold something.

5. Verify the result

For categorical data, the check is exact:

def categories_preserved(src_path, dst_path):
    with rasterio.open(src_path) as s, rasterio.open(dst_path) as d:
        before = set(np.unique(s.read(1, masked=True).compressed()).tolist())
        after = set(np.unique(d.read(1, masked=True).compressed()).tolist())
    invented, lost = sorted(after - before), sorted(before - after)
    print(f"{dst_path:24} invented: {len(invented)} codes {invented[:6]}  lost: {lost or 'none'}")
    return not invented

categories_preserved("landcover.tif", "landcover_bng.tif")
categories_preserved("landcover.tif", "landcover_bilinear.tif")

landcover_bng.tif        invented: 0 codes []  lost: none
landcover_bilinear.tif   invented: 73 codes [11, 12, 13, 14, 15, 16]  lost: none

invented is the failure that matters — a value in the output that was never in the input proves an averaging method was used on labels. lost can be legitimate: a class occupying three pixels may genuinely not survive a tenfold downsample.

For continuous data, compare distributions rather than values:

for name, path in [("source", "elevation.tif"), ("bilinear", "elevation_bng.tif"),
                   ("cubic", "elevation_cubic.tif"), ("no nodata", "elevation_no_nodata.tif")]:
    with rasterio.open(path) as src:
        a = src.read(1, masked=True)
    print(f"{name:9} min {a.min():6.1f}  mean {a.mean():6.2f}  max {a.max():6.1f} m")

source    min   29.0  mean 189.96  max  480.2 m
bilinear  min   29.0  mean 189.98  max  480.1 m
cubic     min   28.3  mean 189.98  max  480.2 m
no nodata min    0.0  mean 188.95  max  480.1 m

The mean barely moves, and the min and max of bilinear stay inside the source's. cubic pushes the minimum outward, from 29.0 to 28.3 m, by overshooting a sharp edge; 0.7 m on a surface with 450 m of relief is why cubic is acceptable on elevation and not on a percentage. The last line is the one to watch for: with no NoData declared, the corners outside the source footprint are written as 0, and a moorland surface acquires a minimum elevation of sea level.

6. Declare NoData, or it spreads

Interpolating methods combine neighbours, so a missing cell has to be known to be missing. Here is the elevation with a 100 × 100 cell hole — a reservoir mask, say — resampled onto a grid shifted by half a cell, three ways:

with rasterio.open("elevation.tif") as src:
    dem, transform, crs = src.read(1), src.transform, src.crs
shifted = transform * Affine.translation(0.5, 0.5)

def shift(data, resampling, **nodata):
    out = np.full(data.shape, np.nan, dtype="float32")
    reproject(data, out, src_transform=transform, src_crs=crs, dst_transform=shifted,
              dst_crs=crs, resampling=resampling, **nodata)
    return out[1:-1, 1:-1]                      # drop the edge that has no source at all

hole = dem.copy()
hole[300:400, 300:400] = np.nan
for label, kwargs in [("declared", {"src_nodata": np.nan, "dst_nodata": np.nan}), ("undeclared", {})]:
    for method in (Resampling.nearest, Resampling.bilinear, Resampling.cubic):
        print(f"NaN {label:10} {method.name:8} {np.isnan(shift(hole, method, **kwargs)).sum():>6,} missing cells")

sentinel = dem.copy()
sentinel[300:400, 300:400] = -9999
out = shift(sentinel, Resampling.bilinear, dst_nodata=np.nan)
blended = (out < -1) & (out > -9999)
print(f"-9999 undeclared, bilinear: {blended.sum()} cells between {out[blended].min():.0f} and {out[blended].max():.0f} m")

NaN declared   nearest  10,000 missing cells
NaN declared   bilinear 10,000 missing cells
NaN declared   cubic    10,000 missing cells
NaN undeclared nearest  10,000 missing cells
NaN undeclared bilinear 10,201 missing cells
NaN undeclared cubic    10,609 missing cells
-9999 undeclared, bilinear: 400 cells between -4929 and -2260 m

With NoData declared, GDAL leaves missing pixels out of every kernel and re-weights the rest, so the hole stays the size it was. Undeclared, a NaN poisons every kernel it touches, and the hole grows with the kernel: 201 more missing cells for bilinear's 2 × 2, 609 for cubic's 4 × 4. Worse is an undeclared sentinel: -9999 is just a very low elevation to the arithmetic, and the cells around the hole come out as plausible-looking but impossible heights thousands of metres below sea level. Always pass src_nodata and dst_nodata.

Changing resolution with out_shape

The data used on this page

Two 0.4° × 0.25° windows over the South Pennines (2.0–1.6°W, 53.6–53.85°N): ESA WorldCover 2021 at 10 m, and the Copernicus GLO-30 elevation model at about 30 m. Both are open, both are read straight from their public cloud buckets as Cloud-Optimised GeoTIFFs, and only the window is downloaded.

import numpy as np
import rasterio
from rasterio.windows import from_bounds

BOUNDS = (-2.0, 53.6, -1.6, 53.85)
SOURCES = {
    "elevation.tif": "https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_W002_00_DEM/"
                     "Copernicus_DSM_COG_10_N53_00_W002_00_DEM.tif",
    "landcover.tif": "https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
                     "ESA_WorldCover_10m_2021_v200_N51W003_Map.tif",
}
CLASSES = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up",
           60: "bare", 70: "snow and ice", 80: "water", 90: "wetland", 95: "mangroves", 100: "moss"}

for name, url in SOURCES.items():
    with rasterio.open(url) as src:
        window = from_bounds(*BOUNDS, transform=src.transform).round_offsets().round_lengths()
        profile = src.profile | {"driver": "GTiff", "width": int(window.width), "height": int(window.height),
                                 "transform": src.window_transform(window), "compress": "deflate",
                                 "tiled": True, "blockxsize": 256, "blockysize": 256}
        with rasterio.open(name, "w", **profile) as dst:
            dst.write(src.read(window=window))
    with rasterio.open(name) as src:
        print(f"{name:14} {src.width} × {src.height}  {src.dtypes[0]:8} nodata {src.nodata}  "
              f"{src.res[0] * 3600:.2f}″ × {src.res[1] * 3600:.2f}″")

elevation.tif  960 × 900  float32  nodata None  1.50″ × 1.00″
landcover.tif  4800 × 3000  uint8    nodata 0.0  0.30″ × 0.30″

The two sources store cells in degrees, and neither kind of cell is square on the ground at 53.7°N — which is why putting them into British National Grid needs resampling at all.

Example 1: changing resolution with out_shape

Rasterio resamples on read when you ask for a different output shape, which avoids a temporary file:

FACTOR = 4        # 30 m → 120 m

with rasterio.open("elevation_bng.tif") as src:
    data = src.read(1, out_shape=(src.height // FACTOR, src.width // FACTOR),
                    resampling=Resampling.average, masked=True)
    # the transform must be scaled to match the new shape
    transform = src.transform * src.transform.scale(src.width / data.shape[-1],
                                                    src.height / data.shape[-2])
    profile = src.profile | {"height": data.shape[0], "width": data.shape[1], "transform": transform}
    forgot = src.profile | {"height": data.shape[0], "width": data.shape[1]}

for name, prof in [("elevation_120m.tif", profile), ("elevation_forgot.tif", forgot)]:
    with rasterio.open(name, "w", **prof) as dst:
        dst.write(data.filled(prof["nodata"]), 1)
    with rasterio.open(name) as check:
        width_km = (check.bounds.right - check.bounds.left) / 1000
        print(f"{name:21} {check.width} × {check.height}  pixel {check.res[0]:.0f} m  covers {width_km:.1f} km east–west")

elevation_120m.tif    220 × 232  pixel 120 m  covers 26.5 km east–west
elevation_forgot.tif  220 × 232  pixel 30 m  covers 6.6 km east–west

src.transform.scale(...) is the part people forget. Reading at a different shape without updating the transform produces a raster that claims 30 m pixels while holding 120 m data — it opens, it plots, and it covers a quarter of the width it should, shrunk into the north-west corner of its real footprint.

Example 2: the same land-cover raster, three ways

WorldCover's codes are spaced ten apart, which makes averaging easy to see. Many land-cover products number their classes 1, 2, 3 … instead:

with rasterio.open("landcover.tif") as src:
    shape = (src.height // 5, src.width // 5)
    results = {m.name: src.read(1, out_shape=shape, resampling=m)
               for m in (Resampling.nearest, Resampling.average, Resampling.mode)}
    raw, nodata, transform = src.read(1), src.nodata, src.transform

for name, arr in results.items():
    print(f"WorldCover codes, {name:8} {len(np.unique(arr)):>3} values: {np.unique(arr)[:10].tolist()}")

dense = {0: 0} | {code: i for i, code in enumerate(sorted(set(np.unique(raw).tolist()) - {0}), start=1)}   # 10 → 1, 20 → 2 …
renumbered = np.vectorize(dense.get)(raw).astype("uint8")
with MemoryFile() as mem:
    with mem.open(driver="GTiff", width=raw.shape[1], height=raw.shape[0], count=1, dtype="uint8",
                  crs="EPSG:4326", transform=transform, nodata=0) as ds:
        ds.write(renumbered, 1)
        avg_dense = ds.read(1, out_shape=shape, resampling=Resampling.average)
print(f"codes 1–{len(dense) - 1},   average  {len(np.unique(avg_dense)):>3} values: {np.unique(avg_dense).tolist()}")

WorldCover codes, nearest    8 values: [10, 20, 30, 40, 50, 60, 80, 90]
WorldCover codes, average   75 values: [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
WorldCover codes, mode       8 values: [10, 20, 30, 40, 50, 60, 80, 90]
codes 1–8,   average    8 values: [1, 2, 3, 4, 5, 6, 7, 8]

With WorldCover's own codes the damage is obvious: average produces dozens of codes that are not classes. With dense codes it produces only codes that exist, because every rounded mean between 1 and 8 is some class's number — the output looks exactly like a valid land-cover map, and a cell that was half tree cover (1) and half built-up (5) comes out as grassland (3). The absence of impossible values does not prove the method was right. Choose by what the data means, not by inspecting the output.

Overviews — the reduced copies a viewer draws when zoomed out — are resampled too, with a method chosen when they are built:

import shutil

for method in (Resampling.nearest, Resampling.average, Resampling.mode):
    shutil.copy("landcover.tif", "landcover_ovr.tif")
    with rasterio.open("landcover_ovr.tif", "r+") as dst:
        dst.build_overviews([2, 4, 8, 16], method)
    with rasterio.open("landcover_ovr.tif") as src:
        zoomed_out = src.read(1, out_shape=(src.height // 16, src.width // 16))   # served from the 1:16 overview
    print(f"overviews built with {method.name:8} {len(np.unique(zoomed_out)):>3} values at 1:16")

overviews built with nearest    8 values at 1:16
overviews built with average   71 values at 1:16
overviews built with mode       6 values at 1:16

A land-cover map with average overviews looks right at full resolution and shows dozens of non-existent classes as soon as you zoom out. mode keeps to real classes but, as in step 3, loses the two rarest; nearest keeps all eight.

Example 3: aligning two rasters onto one grid

Resampling is how you make two rasters comparable cell for cell:

def align_to(src_path, template_path, resampling):
    """Resample src onto the exact grid of template."""
    with rasterio.open(template_path) as tmpl:
        dst_transform, dst_crs, dst_shape = tmpl.transform, tmpl.crs, (tmpl.height, tmpl.width)
    with rasterio.open(src_path) as src:
        nodata = src.nodata if src.nodata is not None else np.nan
        out = np.full(dst_shape, nodata, dtype=src.dtypes[0])
        reproject(rasterio.band(src, 1), out, src_nodata=src.nodata,
                  dst_transform=dst_transform, dst_crs=dst_crs, dst_nodata=nodata,
                  resampling=resampling)
    return out

elevation = align_to("elevation.tif", "landcover_bng.tif", Resampling.bilinear)
with rasterio.open("landcover_bng.tif") as src:
    cover = src.read(1)

print(f"{'class':12} {'cells':>10} {'mean height':>12}")
for code in np.unique(cover[cover != 0]):
    cells = (cover == code) & np.isfinite(elevation)
    print(f"{CLASSES[code]:12} {cells.sum():>10,} {elevation[cells].mean():>10.0f} m")

class             cells  mean height
tree cover    2,049,728        152 m
shrubland           552        368 m
grassland     3,423,886        233 m
cropland        199,208        157 m
built-up      1,615,950        150 m
bare              6,714        166 m
water            36,669        248 m
wetland             660        422 m

Note the two different methods: continuous elevation gets bilinear onto the template, and the categorical template itself was made with nearest in step 4. This is the ordinary case, not an exotic one — aligning a stack almost always needs a per-layer method. The result is also a sanity check on the whole chain, and it matches the ground: towns and woodland sit in the valleys at around 150 m, grassland climbs the hillsides, and water averages higher than grassland because most of the open water here is moorland reservoirs, not rivers. The handful of wetland and shrubland cells are the highest of all, on the tops.

Why grids make resampling unavoidable

Resampling is unavoidable because grids are arbitrary. A raster's cells are an artefact of how it was produced — the sensor's footprint, the processing chain's chosen resolution, the origin someone picked. Two rasters of the same place rarely share a grid, and any operation that combines them, or moves one into a different CRS, must map values from one arrangement of cells to another.

The methods differ in how many source pixels contribute and how their values are combined, and that single distinction explains every behaviour worth knowing.

Nearest neighbour uses one pixel and copies it, so its output value set is a subset of its input value set. That property — not speed, not simplicity — is why it is the correct choice for categorical data. It is a selection, not a computation, and selection is the only operation that respects labels.

Bilinear and cubic compute weighted means, and a weighted mean of category codes is a category code only by accident. This is the same error as taking the mean of postcodes. The reason it survives so long undetected is that the output usually still looks like a plausible raster: it plots, its range is roughly right, its histogram is roughly right, and only a class-by-class comparison reveals that "45" was never a class.

Cubic's overshoot deserves a specific warning. Fitting a smooth curve through four points and evaluating between them can produce values outside the range of those points — the same ringing artefact you see around sharp edges in over-sharpened photographs. On elevation this is harmless. On a bounded quantity — percentage cover, probability, an index constrained to −1…1 — it produces values that are physically impossible, and downstream code that assumed the bound will break or, worse, will not.

Finally, downsampling with a point method is a sampling decision, and so is downsampling with mode. Taking a 10 m raster to 100 m with nearest keeps one pixel in a hundred. That sample is unbiased — step 3 showed every class's share surviving to within a fraction of a percentage point — but noisy, and a small feature survives or vanishes depending on where the output grid happens to fall. mode is the opposite trade: stable cell by cell, systematically biased towards whatever dominates. average for continuous data, and per-class fractions for categorical areas, are the methods that use all the evidence rather than 1% of it.

Upsampling with average or mode, and NoData

• average, mode and sum do nothing special when upsampling. With one source pixel per output pixel there is nothing to aggregate, and the result is identical to nearest.

• NoData only stays put when it is declared. Pass src_nodata and dst_nodata to reproject every time; step 6 shows an undeclared -9999 becoming heights thousands of metres below sea level.

• Cubic and Lanczos on bounded data break the bound. Use cubic_spline if you want smoothness without overshoot, or clip afterwards: np.clip(out, 0, 100).

• Resampling repeatedly compounds error. Reprojecting A→B→C is worse than A→C. Go back to the original raster rather than chaining conversions.

• Overviews have their own resampling method, fixed when they are built (Example 2). rio overview --resampling and gdaladdo -r both default to nearest; the trap is passing average to every file in a folder.

• Resampling.mode is slower than nearest, because it counts. Worth it for a map; see step 3 for why it is the wrong tool for areas.

• Resampling.sum exists for count rasters such as population, where averaging would divide the total.

• gdalwarp -r near|bilinear|cubic|cubicspline|average|mode|sum exposes the same choice on the command line; its default is near.

Internal links

• The raster data model explained — what the values are before you resample them
• How to reproject a raster in Python with Rasterio — the operation that forces this choice most often
• How to merge and mosaic rasters in Python — resampling when grids do not agree
• Raster and vector do not line up in Python — the alignment problem resampling solves
• Rasterio returns the wrong values — NoData interacting badly with interpolation
• How to calculate zonal statistics in Python — where a wrongly resampled grid produces wrong totals
• Introduction to Rasterio — reading and writing basics
• Coordinate precision and floating point in GIS explained — why grid alignment is so fragile
• Fragmentation metrics change with resolution — what aggregating a forest map does to patches

FAQ

Which resampling method should I use by default?

There is no safe default. For continuous data use bilinear (or average when downsampling); for categorical data use nearest (or mode when downsampling for a map, and class fractions when you need areas). If you do not know which you have, nearest is the only method that cannot invent values.

How do I know whether my raster is categorical?

Count distinct values. A handful of integers in an integer dtype means labels. Hundreds of thousands of values, or any float dtype, means measurements. Elevation in whole metres is continuous despite being an integer — count first.

Why did my land-cover classes change after reprojection?

An interpolating method averaged the class codes. Reproject again with Resampling.nearest from the original file; do not try to repair the corrupted output.

Is cubic better than bilinear?

Smoother, a little slower, and able to produce values outside the input range near sharp edges — on the tree-cover percentage in step 2, 14% of cells landed outside 0–100. Fine for elevation and hillshades; for anything with a hard physical bound use bilinear or cubic_spline, which does not overshoot.

What method should overviews use?

The same logic as the data: average for continuous, mode or nearest for categorical. rio overview --build 2^1..4 --resampling average landcover.tif produces a zoomed-out view that shows classes which do not exist.

Does resampling lose data?

Downsampling always does. Upsampling adds no information and can add artefacts. Resampling twice compounds both, so always resample from the original.

Why do my two aligned rasters still disagree by half a pixel?

The output grid origin was probably computed rather than taken from a template. Reproject both onto the same explicit transform, as in Example 3, rather than letting each one pick its own default grid.