This guide is part of Raster Alignment & Resampling Techniques, itself a section of the wider Automated Vector & Raster Cleaning Workflows reference.

Choosing Resampling Methods for Categorical and Continuous Rasters

The Problem: The Wrong Kernel Silently Fabricates Data

Every resample either interpolates between neighbouring pixels or copies one of them, and applying an interpolating kernel to categorical data invents class codes that never existed in the source.

This is one of the most damaging silent errors in raster ETL, because the output still looks like valid data:

  • Fabricated class codes — bilinear on a land-cover raster averages class 3 and class 5 into a 4, which may mean an unrelated class; every subsequent area tally and reclassification inherits the error.
  • Broken class boundaries — interpolation smears sharp edges between land-cover types into a fringe of invalid intermediate values along every boundary.
  • Corrupted joins and zonal statistics — a masked training raster resampled with the wrong kernel no longer matches its imagery, poisoning supervised classification and per-class statistics.
  • Aliasing on continuous downscale — the opposite mistake, using nearest to shrink a DEM, drops entire rows and columns and produces a jagged, aliased surface instead of a smoothed one.

The remedy is a small decision function that maps data type and direction of scaling to a safe kernel — and refuses to interpolate categorical data at all.

Version and Environment Compatibility

rasterio GDAL backend Kernels available for reads Known caveat
≥ 1.3.9 ≥ 3.6 nearest, bilinear, cubic, cubic_spline, lanczos, average, mode Stable; all read-resampling kernels honoured
1.3.0–1.3.8 3.4–3.5 same set mode slow on very large decimation windows
1.2.x 3.2–3.3 nearest, bilinear, cubic, average, mode Resampling.sum/rms unavailable for reads
< 1.2 < 3.2 Not recommended Limited read-resampling support; upgrade required 
pip install "rasterio>=1.3.9" "numpy>=1.26.0"

Decision Logic: Data Type to Resampling Kernel

Resampling kernel decision tree by data type and scale direction A decision tree starting from the raster data type. Categorical rasters branch to mode when downscaling and nearest when upscaling. Continuous rasters branch to average when downscaling and bilinear or cubic when upscaling. Raster data type? dtype + intent Categorical land cover · masks Continuous elevation · reflectance mode downscale majority class nearest upscale copies class average downscale area-weighted mean bilinear / cubic upscale smooth gradient

The table below is the same logic in reference form. The direction of scaling matters: mode and average aggregate many source pixels into one when downscaling, while nearest, bilinear, and cubic sample from few pixels when upscaling.

Data type Example Downscale (scale < 1) Upscale (scale > 1) Never use
Categorical Land cover, classified mask Resampling.mode Resampling.nearest bilinear, cubic, lanczos, average
Continuous (float) Elevation, reflectance, temperature Resampling.average Resampling.bilinear / cubic mode
Continuous (scaled int) Int16 scaled reflectance Resampling.average Resampling.bilinear mode, nearest
Count / density Population per pixel Resampling.sum Resampling.bilinear mode

Recipe: resample_raster with a Categorical Guard

The helper maps data type and scale factor to a Resampling kernel, refuses interpolating kernels on categorical rasters, and resamples with a scaled read. The categorical flag lets callers override the integer-means-categorical default for scaled-integer continuous data.

import logging
from pathlib import Path

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

logger = logging.getLogger(__name__)

# Kernels that interpolate between neighbouring pixels and therefore fabricate
# values that are not present in the source — unsafe on categorical data.
_INTERPOLATING: set[Resampling] = {
    Resampling.bilinear,
    Resampling.cubic,
    Resampling.cubic_spline,
    Resampling.lanczos,
    Resampling.average,
}


def recommend_resampling(dtype: str, categorical: bool, scale: float) -> Resampling:
    """Map a data type and scale direction to a safe Resampling kernel."""
    if categorical:
        # Aggregate to the majority class when shrinking; copy when enlarging.
        return Resampling.mode if scale < 1.0 else Resampling.nearest
    # Continuous: area-weighted mean avoids aliasing on downscale.
    return Resampling.average if scale < 1.0 else Resampling.bilinear


def resample_raster(
    src_path: Path,
    out_path: Path,
    scale: float,
    method: Resampling | None = None,
    categorical: bool | None = None,
) -> Resampling:
    """
    Resample a raster by a scale factor, choosing or validating the kernel
    against the raster's data type.

    Parameters
    ----------
    src_path:    Source raster on disk.
    out_path:    Destination GeoTIFF path.
    scale:       Output size relative to input; 0.5 halves resolution,
                 2.0 doubles it. Must be positive.
    method:      Explicit Resampling kernel. If None, one is chosen from the
                 data type and scale direction.
    categorical: Override for classification. Defaults to True for integer
                 dtypes and False for floating-point dtypes.

    Returns
    -------
    The Resampling kernel actually applied.
    """
    if scale <= 0:
        raise ValueError(f"scale must be positive, got {scale}")

    with rasterio.open(src_path) as src:
        dtype = src.dtypes[0]
        is_categorical = (
            categorical
            if categorical is not None
            else np.issubdtype(np.dtype(dtype), np.integer)
        )

        if method is None:
            method = recommend_resampling(dtype, is_categorical, scale)
            logger.info(
                "Auto-selected %s (dtype=%s, categorical=%s, scale=%.3f)",
                method, dtype, is_categorical, scale,
            )

        # Guard: interpolating a categorical raster invents invalid class codes.
        if is_categorical and method in _INTERPOLATING:
            raise ValueError(
                f"{method} interpolates and will fabricate invalid class codes on "
                f"categorical raster {src_path}. Use Resampling.nearest or "
                f"Resampling.mode instead."
            )

        new_height = max(1, int(round(src.height * scale)))
        new_width = max(1, int(round(src.width * scale)))

        # Read into a scaled out_shape; GDAL applies the kernel during the read.
        data = src.read(
            out_shape=(src.count, new_height, new_width),
            resampling=method,
        )

        # Scale the affine so real-world coordinates stay fixed.
        dst_transform = src.transform * src.transform.scale(
            src.width / new_width,
            src.height / new_height,
        )

        profile = src.profile.copy()
        profile.update(
            height=new_height,
            width=new_width,
            transform=dst_transform,
            compress="deflate",
            tiled=True,
            blockxsize=512,
            blockysize=512,
        )

        out_path.parent.mkdir(parents=True, exist_ok=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(data)

    logger.info(
        "Resampled %s by scale=%.3f -> %dx%d with %s",
        src_path, scale, new_width, new_height, method,
    )
    return method

Key Implementation Notes

  • Integer means categorical by default, but the caller can override. Most uint8/int16 rasters are classes or masks, so defaulting them to categorical is the safe assumption. Scaled reflectance and quantised elevation are the exceptions — pass categorical=False so they get average or bilinear rather than mode.

  • The guard fires before any pixel is read. Raising on an interpolating kernel plus a categorical raster up front turns a silent data-corruption bug into an immediate, obvious failure at the call site instead of a wrong map three stages later.

  • average is treated as interpolating for the guard. Even though it aggregates rather than fits a curve, average still produces fractional values that are not valid class codes, so it belongs in the blocked set for categorical data alongside bilinear and cubic.

  • Scale direction, not just data type, drives the recommendation. mode and average only make sense when shrinking, where each output pixel summarises many inputs. When enlarging, nearest copies the class and bilinear smooths the surface.

  • The transform is rescaled with transform.scale. Resampling changes pixel size, so the affine must change with it or the raster’s real-world footprint shifts. Multiplying by transform.scale(width_ratio, height_ratio) keeps the extent fixed while the pixel count changes.

  • Resampling.mode can be slow on heavy decimation. Reducing a large raster by a big factor forces mode to tally every source pixel in each output window. For extreme downscaling of categorical data, build overviews once and read from them instead of resampling in a hot loop.

Troubleshooting Kernel Selection Failures

Failure Root Cause Fix
Land-cover raster gains classes that never existed Bilinear or cubic applied to categorical data Use Resampling.nearest or Resampling.mode; rely on the guard
ValueError: ... fabricate invalid class codes Guard caught an interpolating kernel on categorical data Pass a categorical-safe kernel, or set categorical=False if the ints are continuous
Downscaled DEM looks jagged and aliased nearest used to shrink a continuous surface Use Resampling.average for downscaling continuous rasters
Scaled-integer reflectance loses detail as blocky classes Integer dtype auto-classified as categorical Pass categorical=False so continuous kernels apply
Resampled raster shifts off its footprint Transform not rescaled to the new pixel size Multiply by src.transform.scale(width_ratio, height_ratio)

Integration Note

Kernel selection is the decision that sits inside the warp step of Raster Alignment & Resampling Techniques: whenever a raster changes resolution or projection, recommend_resampling should choose the kernel before the pixels move. When you push several rasters onto one shared grid with the reproject_to_grid recipe for a common grid, pass the kernel this helper returns as that function’s resampling argument so categorical and continuous inputs are each warped correctly onto the same lattice.