Bilinear vs Nearest-Neighbour Resampling for Ecological Covariates

Every time you reproject or regrid an ecological covariate raster, GDAL has to invent values at the new pixel centres, and the interpolation rule you pass decides whether those values are physically meaningful or quietly corrupted. Pick Resampling.bilinear for a categorical land-cover grid and you manufacture class codes that never existed; pick Resampling.nearest for a continuous elevation surface and you inject blocky terracing that propagates into every slope, aspect, and index derived from it. This guide narrows in on that single choice — bilinear versus nearest-neighbour (and when cubic earns its keep) — for the covariate stacks that feed forest and habitat models. It sits inside the Coordinate Reference Systems for Forestry workflow, part of the broader Ecological GIS Data Foundations in Python framework for reproducible spatial analysis.

When to Use Each Resampling Method

The decision is driven almost entirely by one property of the raster: is the cell value a measurement on a continuous scale (temperature, elevation, reflectance, canopy height) or a label from a finite set of classes (land-cover type, soil group, sampling stratum)? Interpolating between two class codes — say, 41 (deciduous forest) and 52 (shrub) — to produce 46 is nonsense; there is no half-forest-half-shrub category. Averaging two elevations, by contrast, is exactly what you want.

The two rows that cover the vast majority of covariate work are the first two. A rule of thumb: if you would ever call numpy.unique() on the raster and expect a short list of integer codes, use nearest; otherwise use bilinear. Cubic convolution is worth the extra cost only for continuous surfaces where you care about crisp gradients — a digital elevation model destined for terrain derivatives, or a canopy height model — and where a small out-of-range overshoot is acceptable or clamped afterwards.

Downsampling (going to a coarser pixel) is a special case: nearest and bilinear both throw away most of the contributing cells, so prefer mode for categorical and average for continuous when the target resolution is markedly coarser than the source. These continuous resampling choices flow directly into the predictor grids assembled in Environmental Predictor Stacking, where every layer must share one grid before modelling.

Minimal Reproducible Example

The snippet reprojects two covariates from their native CRS to a common analytical grid (EPSG:5070, CONUS Albers equal-area): a continuous mean-temperature raster with bilinear, and a categorical land-cover raster with nearest. Both use rasterio.warp.reproject driven by the Resampling enum, and a single helper enforces the method per raster type.

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

DST_CRS = "EPSG:5070"  # NAD83 / CONUS Albers — equal-area, metres


def warp_raster(src_path: str, dst_path: str, resampling: Resampling) -> None:
    """Reproject a single-band raster to DST_CRS with an explicit resampling rule."""
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, DST_CRS, src.width, src.height, *src.bounds
        )
        meta = src.meta.copy()
        meta.update(crs=DST_CRS, transform=transform, width=width, height=height)

        with rasterio.open(dst_path, "w", **meta) as dst:
            reproject(
                source=rasterio.band(src, 1),
                destination=rasterio.band(dst, 1),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=DST_CRS,
                resampling=resampling,
                src_nodata=src.nodata,
                dst_nodata=src.nodata,
            )


# Continuous covariate: smooth interpolation is correct.
warp_raster("mean_temp.tif", "mean_temp_5070.tif", Resampling.bilinear)

# Categorical covariate: never interpolate class codes.
warp_raster("landcover.tif", "landcover_5070.tif", Resampling.nearest)

The only thing that changes between the two calls is the Resampling member. That is the whole point: the reprojection machinery is identical, and the correctness of the output hinges on matching the interpolation rule to the data semantics. For a stack of many layers, drive it from a small manifest so the method is recorded, not guessed:

COVARIATES = [
    ("mean_temp.tif",   "mean_temp_5070.tif",   Resampling.bilinear),
    ("elevation.tif",   "elevation_5070.tif",   Resampling.cubic),
    ("landcover.tif",   "landcover_5070.tif",   Resampling.nearest),
    ("soil_class.tif",  "soil_class_5070.tif",  Resampling.nearest),
]

for src_path, dst_path, method in COVARIATES:
    warp_raster(src_path, dst_path, method)

Parameter Reference

These are the reproject arguments that change the ecological correctness of the output, not just its mechanics.

A subtlety that bites forestry pipelines: if src_nodata is unset on a continuous raster, bilinear will happily average the nodata sentinel (often -9999) with valid neighbours along every masked edge, dragging border cells toward large negative values. Always declare nodata before warping.

Expected Output and Verification

The defining, checkable property of a correct categorical resample is that no new class codes appear: the set of output values must be a subset of the input values (plus nodata). For a continuous resample, the output range must stay within the input range for nearest and bilinear (cubic may slightly overshoot), and the smoothing must not shift the mean appreciably.

import numpy as np
import rasterio


def assert_categorical_preserved(src_path: str, dst_path: str) -> None:
    """A nearest resample must invent no new class codes."""
    with rasterio.open(src_path) as s:
        src_vals = set(np.unique(s.read(1)).tolist())
    with rasterio.open(dst_path) as d:
        dst_vals = set(np.unique(d.read(1)).tolist())
    new_codes = dst_vals - src_vals
    assert not new_codes, f"Resampling invented class codes: {sorted(new_codes)}"


def assert_continuous_bounded(src_path: str, dst_path: str, tol: float = 1e-3) -> None:
    """A bilinear resample must not exceed the source value range."""
    with rasterio.open(src_path) as s:
        a = s.read(1, masked=True)
    with rasterio.open(dst_path) as d:
        b = d.read(1, masked=True)
    assert b.min() >= a.min() - tol, "Output dipped below source minimum"
    assert b.max() <= a.max() + tol, "Output exceeded source maximum — used cubic?"


assert_categorical_preserved("landcover.tif", "landcover_5070.tif")
assert_continuous_bounded("mean_temp.tif", "mean_temp_5070.tif")

If assert_categorical_preserved fails, you almost certainly warped a class raster with bilinear. If a continuous check fails only on the maximum, you likely used cubic, whose convolution kernel can overshoot at sharp edges — clamp with np.clip to the physical range if that matters (for example, a reflectance layer bounded to [0, 1]). The verified continuous grids are then ready to feed a Vegetation Index Calculation in Python step, where any resampling artefact would propagate straight into NDVI or EVI.

Common Pitfalls

• Interpolating class codes. Warping land-cover, soil, or strata rasters with bilinear/cubic blends integer codes into fractional values that map to the wrong — or nonexistent — class. Categorical always means nearest (upsampling) or mode (downsampling).
• Terracing a DEM with nearest. Reprojecting continuous elevation or canopy height with nearest copies the same value across blocks, creating stair-step artefacts that inflate slope and aspect at every block edge. Use bilinear or cubic for terrain.
• Nodata leaking into the average. With src_nodata unset, bilinear averages the fill sentinel into edge cells, producing a fringe of physically impossible values. Set src_nodata and dst_nodata explicitly.
• Resampling then stacking mismatched grids. Warping each covariate to the same CRS but not the same transform and resolution leaves layers one pixel offset; pin dst_resolution and align the transform so the stack is cell-for-cell coincident.

Frequently Asked Questions

Related

• Coordinate Reference Systems for Forestry — the parent workflow where reprojection and resampling choices are set.
• Vegetation Index Calculation in Python — spectral indices whose inputs must be resampled without artefacts.
• Environmental Predictor Stacking — assembling many resampled covariates onto one shared grid.
• Ecological GIS Data Foundations in Python — the reproducible-pipeline framework this guide belongs to.