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

Problem statement #

You reproject a land-cover raster and the class counts change:

before = np.unique(src.read(1), return_counts=True)
# classes: [1 2 3 4 5]  counts: [ 4021  9330  1188 20044  5417]

after = np.unique(dst.read(1), return_counts=True)
# classes: [1 2 3 4 5]  counts: [ 4106  9241  1301 19832  5520]

Slightly different, which is expected. Then you look closer:

print(np.unique(dst.read(1)))

[1 2 3 4 5]

Fine here. But run the same reprojection with the default settings on a raster whose classes are 10, 20, 30, 40 and you get:

[10 12 14 17 20 23 25 28 30 34 37 40]

Classes that do not exist. Nothing errored. The resampling method averaged category codes, and "class 17" is the arithmetic mean of "deciduous forest" and "water".

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, and choosing the default is choosing wrong roughly half the time.

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,
    dst_transform=dst_transform, dst_crs="EPSG:27700",
    resampling=Resampling.nearest,     # ← never leave this to chance
)

Step-by-step solution #

1. Work out what your numbers mean #

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

import numpy as np, rasterio

with rasterio.open("layer.tif") as src:
    data = src.read(1, masked=True)

vals = np.unique(data.compressed())
print(f"{len(vals)} distinct values, dtype {src.dtypes[0]}")
print(vals[:12])

7 distinct values, dtype uint8
[1 2 3 4 5 6 7]

Seven distinct integers in a uint8 raster is categorical, near-certainly. Compare:

418293 distinct values, dtype float32
[-3.21  -3.19  -3.18 ... ]

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

Resampling.nearest

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

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

Resampling.bilinear

• 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 cubic spline use a 4 × 4 neighbourhood.

Resampling.cubic

• Smoother still, better at preserving gradients.
• Can overshoot: near a sharp step, output values can fall outside the input range. On a 0–100 percentage raster, cubic can produce 103.
• About four times the cost of bilinear.

Average and mode are aggregating methods, only meaningful when several source pixels fall inside one output pixel.

Resampling.average   # mean of contributing pixels — continuous data
Resampling.mode      # most common contributing value — categorical data

• These are the right answers for downsampling, and both are ignored when upsampling.
• average on a masked raster averages only valid pixels in recent GDAL versions; check yours if edge cells matter.

3. Match the method to the direction #

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

Downsampling (coarser output) means many input pixels per output pixel. Nearest neighbour throws away all but one of them, so a 10 m raster taken to 100 m keeps 1% of the information and discards 99% — including, quite possibly, every pixel of the feature you cared about. Use average for continuous data and mode for categorical.

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.

4. Set it explicitly, every time #

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

def reproject_raster(src_path, dst_path, dst_crs, resampling):
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds
        )
        profile = src.profile.copy()
        profile.update(crs=dst_crs, transform=transform, width=width, height=height)

        with rasterio.open(dst_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform, src_crs=src.crs,
                    dst_transform=transform, dst_crs=dst_crs,
                    resampling=resampling,
                )

reproject_raster("landcover.tif", "landcover_bng.tif",
                 "EPSG:27700", Resampling.nearest)
reproject_raster("elevation.tif", "elevation_bng.tif",
                 "EPSG:27700", Resampling.bilinear)

Two calls, two methods, and the difference is not cosmetic.

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()))
        after = set(np.unique(d.read(1, masked=True).compressed()))
    invented = after - before
    lost = before - after
    print(f"invented: {sorted(invented) or 'none'}")
    print(f"lost:     {sorted(lost) or 'none'}")
    return not invented

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

invented: none
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 [("before", "elevation.tif"), ("after", "elevation_bng.tif")]:
    with rasterio.open(path) as src:
        a = src.read(1, masked=True)
    print(f"{name:7} min {a.min():7.1f}  mean {a.mean():7.2f}  max {a.max():7.1f}")

before     -3.0   214.83   1344.0
after      -3.0   214.79   1344.0

The mean should barely move. The min and max should not move outward — if they do, cubic overshoot has invented values beyond the real range.

Code examples #

Example 1: changing resolution with out_shape #

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

import rasterio
from rasterio.enums import Resampling

FACTOR = 4        # 10 m → 40 m

with rasterio.open("elevation_10m.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.copy()

profile.update(height=data.shape[0], width=data.shape[1], transform=transform)

with rasterio.open("elevation_40m.tif", "w", **profile) as dst:
    dst.write(data.filled(profile["nodata"]), 1)

src.transform.scale(...) is the part people forget. Reading at a different shape without updating the transform produces a raster that claims 10 m pixels while holding 40 m data — it opens, it plots, and it is wrong by a factor of four.

Example 2: the same land-cover raster, both ways #

import numpy as np, rasterio
from rasterio.enums import Resampling

with rasterio.open("landcover_10m.tif") as src:
    shape = (src.height // 5, src.width // 5)
    near = src.read(1, out_shape=shape, resampling=Resampling.nearest)
    avg  = src.read(1, out_shape=shape, resampling=Resampling.average)
    mode = src.read(1, out_shape=shape, resampling=Resampling.mode)

for name, arr in [("nearest", near), ("average", avg), ("mode", mode)]:
    print(f"{name:8} {sorted(np.unique(arr))}")

nearest  [1, 2, 3, 4, 5]
average  [1, 2, 3, 4, 5]
mode     [1, 2, 3, 4, 5]

With class codes 1–5 the damage from average is invisible in the class list, because rounding to uint8 lands back on real codes. That is what makes this bug so persistent. Now with realistic codes:

lookup = {1: 10, 2: 20, 3: 30, 4: 40, 5: 50}
coded = np.vectorize(lookup.get)(near).astype("uint8")
# ... resample `coded` with average ...
print(sorted(np.unique(avg_coded)))

[10, 13, 16, 20, 24, 27, 30, 33, 37, 40, 44, 50]

Now the corruption is obvious. The absence of impossible values does not prove the method was right — it can just mean the class codes are dense enough to hide it. Choose by what the data means, not by inspecting the output.

Example 3: aligning two rasters onto one grid #

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

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

def align_to(src_path, template_path, resampling):
    """Resample src onto the exact grid of template."""
    with rasterio.open(template_path) as tmpl:
        dst_profile = tmpl.profile.copy()
        dst_shape = (tmpl.height, tmpl.width)
        dst_transform, dst_crs = tmpl.transform, tmpl.crs

    with rasterio.open(src_path) as src:
        dst_profile.update(dtype=src.dtypes[0], nodata=src.nodata, count=src.count)
        out = np.empty((src.count, *dst_shape), dtype=src.dtypes[0])
        reproject(
            source=src.read(),
            destination=out,
            src_transform=src.transform, src_crs=src.crs,
            src_nodata=src.nodata,
            dst_transform=dst_transform, dst_crs=dst_crs,
            dst_nodata=src.nodata,
            resampling=resampling,
        )
    return out, dst_profile

slope, profile = align_to("slope_30m.tif", "landcover_10m.tif", Resampling.bilinear)
cover, _       = align_to("landcover_25m.tif", "landcover_10m.tif", Resampling.nearest)

# now these two arrays are cell-for-cell comparable
mean_slope_by_class = {
    int(c): float(slope[0][cover[0] == c].mean()) for c in np.unique(cover[0])
}

Note the two different methods in the two calls — continuous slope gets bilinear, categorical cover gets nearest. This is the ordinary case, not an exotic one: aligning a stack almost always needs a per-layer method.

Explanation #

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 "17" 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 error, not an approximation. Taking a 10 m raster to 100 m with nearest keeps one pixel in a hundred. If the feature of interest occupies a small share of the area — buildings, watercourses, roads — the output systematically under-represents it, and by an amount that depends on where the output grid happens to fall. average for continuous data and mode for categorical are not merely nicer; they are the methods that use all the evidence rather than 1% of it.

Edge cases or notes #

• Resampling.average and Resampling.mode are ignored when upsampling. With one source pixel per output pixel there is nothing to aggregate, and the result matches nearest.
• NoData spreads under interpolation. A bilinear output cell touching one NoData source cell becomes NoData (or contaminated, depending on GDAL version). Pass src_nodata and dst_nodata explicitly to reproject.
• Cubic on bounded data can break the bound. 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 use their own resampling method. rio overview --resampling average — the default for a categorical raster produces a wrong-looking zoomed-out view.
• Resampling.mode is slow relative to nearest, because it counts. On a large downsample this is worth it; on a small one it is not.
• A "sum" method exists (Resampling.sum, GDAL 3.1+) for count rasters such as population, where averaging would destroy the total.
• gdalwarp -r near|bilinear|cubic|average|mode exposes the same choice on the command line, with the same defaults problem.

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

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). 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, four times slower, and able to produce values outside the input range near sharp edges. Better for visual products like hillshades; not better for anything with a hard physical bound.

What method should overviews use? #

The same logic as the data: average for continuous, mode or nearest for categorical. rio overview --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.