Resampling with GDAL Warp vs rasterio reproject

Both call the same warp engine; the results match when the options match:

gdalwarp -t_srs EPSG:32636 -tr 10 10 -tap -r bilinear \
         -srcnodata 0 -dstnodata 0 -multi -wo NUM_THREADS=ALL_CPUS \
         -co TILED=YES -co COMPRESS=ZSTD in.tif out.tif
reproject(rasterio.band(src, 1), rasterio.band(dst, 1),
          dst_crs="EPSG:32636", dst_transform=aligned_transform,
          resampling=Resampling.bilinear, src_nodata=0, dst_nodata=0,
          num_threads=os.cpu_count())

The choice is about where the step lives in a pipeline, not about output quality. This page belongs to handling pixel resolution and scaling in Core Raster Fundamentals & STAC Mapping.

One Engine, Three Front Ends

Same warp engine underneath The gdalwarp command-line tool, rasterio's warp.reproject function and rioxarray's rio.reproject method are three interfaces to one GDAL warp engine. They differ in defaults and in how options are expressed, which is where output differences come from, not in the resampling mathematics. Three interfaces, one implementation gdalwarp shell, file to file rasterio reproject arrays or band handles rioxarray reproject labelled DataArrays GDAL warp engine same kernels, same approximations, same threads Differences between outputs are differences in options, never in mathematics.

Knowing that removes a surprising amount of confusion. When two tools produce different outputs from the same input, the question is never “which algorithm is better” but “which option did one of them set that the other did not” — and the answer is almost always one of three: target alignment, the approximation threshold, or nodata.

Environment & Setup

Package Version pin Used for
GDAL >=3.6 gdalwarp and the shared warp engine
rasterio >=1.3.0 warp.reproject and calculate_default_transform
rioxarray >=0.15 rio.reproject on labelled arrays
numpy >=1.23 Comparing outputs 
pip install "rasterio>=1.3.0" "rioxarray>=0.15" "numpy>=1.23"

Complete Working Example

import os
import subprocess

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject, transform_bounds


def aligned_grid(src_path: str, dst_crs: str, res: float):
    """Target transform and shape snapped to multiples of res — gdalwarp's -tap."""
    with rasterio.open(src_path) as src:
        l, b, r, t = transform_bounds(src.crs, dst_crs, *src.bounds, densify_pts=21)
    l, b = np.floor(l / res) * res, np.floor(b / res) * res
    r, t = np.ceil(r / res) * res, np.ceil(t / res) * res
    return from_origin(l, t, res, res), int(round((t - b) / res)), int(round((r - l) / res))


def warp_python(src_path: str, dst_path: str, dst_crs: str, res: float,
                method: Resampling = Resampling.bilinear) -> None:
    transform, h, w = aligned_grid(src_path, dst_crs, res)
    with rasterio.open(src_path) as src:
        profile = src.profile | {"crs": dst_crs, "transform": transform, "width": w,
                                 "height": h, "tiled": True, "compress": "zstd"}
        with rasterio.open(dst_path, "w", **profile) as dst:
            for b in range(1, src.count + 1):
                reproject(rasterio.band(src, b), rasterio.band(dst, b),
                          src_crs=src.crs, src_transform=src.transform,
                          dst_crs=dst_crs, dst_transform=transform,
                          src_nodata=src.nodata, dst_nodata=src.nodata,
                          resampling=method, num_threads=os.cpu_count(),
                          warp_mem_limit=512, init_dest_nodata=True)


def warp_cli(src_path: str, dst_path: str, dst_crs: str, res: float,
             method: str = "bilinear") -> None:
    subprocess.run([
        "gdalwarp", "-overwrite", "-t_srs", dst_crs, "-tr", str(res), str(res), "-tap",
        "-r", method, "-multi", "-wo", "NUM_THREADS=ALL_CPUS", "-wm", "512",
        "-co", "TILED=YES", "-co", "COMPRESS=ZSTD", src_path, dst_path,
    ], check=True)


if __name__ == "__main__":
    warp_python("B11_20m.tif", "B11_py.tif", "EPSG:32636", 10.0)
    warp_cli("B11_20m.tif", "B11_cli.tif", "EPSG:32636", 10.0)
    with rasterio.open("B11_py.tif") as a, rasterio.open("B11_cli.tif") as b:
        assert a.transform == b.transform, (a.transform, b.transform)
        diff = np.abs(a.read(1).astype("float64") - b.read(1).astype("float64"))
        print("max abs difference:", diff.max())

aligned_grid reproduces what -tap does in the command-line tool: it snaps the target bounds to multiples of the resolution so that outputs from different inputs share a grid. Omitting it in the Python path while using -tap on the command line is the single most common reason the two outputs differ — by a fraction of a pixel in origin, which makes every value differ slightly.

Where Outputs Diverge

Defaults that differ and what they change Target alignment differs because gdalwarp aligns only with -tap and rasterio uses whatever transform it is given. The approximation error threshold defaults to an eighth of a pixel in gdalwarp and can differ in library calls, changing values slightly along curved warps. Nodata handling differs when one side passes source and destination nodata and the other does not, changing edges. Three settings explain almost every mismatch option symptom when mismatched target alignment (-tap) every pixel differs a little; transforms differ in origin error threshold (-et) tiny differences growing toward the edges of a large warp source/destination nodata a one-pixel band of differences along every nodata edge Match these three and the outputs agree to floating-point precision.

The error threshold is the subtlest. GDAL approximates the coordinate transformation with linear interpolation across small regions and recomputes exactly only when the approximation error exceeds a threshold — an eighth of a pixel by default in gdalwarp. That is invisible in almost every use, but it means two runs with different thresholds give values differing in the last decimal places, especially across large, strongly curved warps. Set it explicitly on both sides if byte-identical output matters, or to zero for an exact transform at some cost in speed.

Choosing Between Them

The tools are interchangeable in output, so the choice is about pipeline shape. Use the command-line tool when the step is a file-to-file conversion with nothing downstream in Python — a batch job converting an archive, a shell script, a container entry point. Use rasterio.warp.reproject when the result is consumed as arrays immediately afterwards, because warping directly into a NumPy destination avoids writing and re-reading an intermediate file. Use rioxarray when the data is already a labelled DataArray and you want the coordinates and CRS carried automatically, as in aligning two rasters with reproject_match.

Whichever you pick, express the target grid explicitly rather than letting the tool derive it. A grid derived from each input’s own bounds differs from file to file, so outputs that should stack do not — the underlying reason -tap and explicit transforms exist at all.

Verification

Transforms equal, values within tolerance Two warp outputs are compared first by transform, which must be identical, and then by maximum absolute value difference, which should be zero or at floating-point noise level when options are matched. A matching transform with differences concentrated at nodata edges points to nodata handling; differences everywhere point to alignment. Diagnose a mismatch in two steps 1 · transforms equal? no → alignment option differs between the runs 2 · where are the diffs? at nodata edges → nodata growing outward → threshold Map the difference raster; its pattern names the unmatched option. 
import numpy as np
import rasterio

with rasterio.open("B11_py.tif") as a, rasterio.open("B11_cli.tif") as b:
    assert a.transform == b.transform and a.shape == b.shape
    da, db = a.read(1, masked=True), b.read(1, masked=True)
    assert (da.mask == db.mask).all(), "nodata footprints differ"
    diff = np.abs(da.astype("float64") - db.astype("float64"))
    print(f"max {diff.max():.6f}, mean {diff.mean():.8f}")

Mapping the difference raster rather than only summarising it is the fastest diagnosis: a stripe along nodata edges names the nodata option, a gradient growing outward names the error threshold, and a uniform speckle across the whole image names alignment.

Common Errors

Outputs differ everywhere by small amounts

The target grids are offset by a fraction of a pixel. Use -tap on the command line and an explicitly snapped transform in Python.

A one-pixel halo of differences along edges

One side passed nodata to the warp and the other did not. Pass source and destination nodata to both.

The Python warp is much slower

num_threads was left at one. Pass it explicitly, along with a warp_mem_limit large enough for the chunk.

Classes change after resampling

Bilinear or cubic was used on categorical data. Use nearest or mode, as covered in choosing the right resampling method for Sentinel-2.

Frequently Asked Questions

Q: Do gdalwarp and rasterio reproject give the same result? Yes, when every option is matched — both call the same GDAL warp engine. Differences almost always trace to one side using a default the other overrides: target alignment, the approximation error threshold, or nodata handling.

Q: Which is faster? They are the same engine, so equally fast with equal settings. gdalwarp exposes multithreading and warp memory flags directly; in rasterio pass num_threads and warp_mem_limit to reproject to get the same behaviour.

Q: When is the command-line tool the better choice? For one-off conversions, shell pipelines and batch jobs where no further Python processing follows. For anything that reads the result back into arrays, rasterio avoids a round trip through disk and keeps the whole step testable in Python.

Q: Can GDAL’s Python bindings be used instead? Yes — gdal.Warp accepts the same options as the command line as keyword arguments and runs in-process. It is a reasonable middle ground when you want gdalwarp’s option names inside Python without shelling out.