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

Reprojecting Rasters to a Common Grid with rasterio.warp

The Problem: Independently Reprojected Rasters Never Truly Line Up

Reprojecting each raster on its own — even to the same CRS and the same nominal resolution — does not guarantee that the results share a grid, because each output’s pixel origin is rounded independently from its own bounding box.

This breaks raster stacking and band math in ways that are easy to miss:

  • ValueError on np.stack — two outputs that differ by a single row or column cannot be stacked into one array, halting the pipeline at the combine step rather than at ingestion.
  • Sub-pixel misregistration in band math — when origins differ by a fraction of a pixel, a difference or ratio between two “aligned” bands compares slightly different ground locations, producing plausible-looking but wrong results.
  • Non-integer affine coefficients — origins that do not fall on a clean lattice break tile-aligned reads of Cloud-Optimized GeoTIFFs and force GDAL to re-block on every read.
  • Irreproducible outputs — floating-point drift in per-source calculate_default_transform calls makes identical inputs yield subtly different grids between runs, defeating content hashing and idempotency checks.

The fix is to treat one target CRS, resolution, and bounds as a fixed contract, derive the destination transform deterministically from that contract, and reproject every source into it.

Version and Environment Compatibility

rasterio GDAL backend Recommended API Known caveat
≥ 1.3.9 ≥ 3.6 reproject + transform.from_origin Stable; num_threads honoured for warp
1.3.0–1.3.8 3.4–3.5 reproject + from_origin warp_mem_limit may be ignored on Windows
1.2.x 3.2–3.3 reproject + from_origin calculate_default_transform ignores dst_bounds; build the grid from bounds manually
< 1.2 < 3.2 Not recommended Resampling.sum/rms unavailable; upgrade required 
pip install "rasterio>=1.3.9" "numpy>=1.26.0" "affine>=2.4.0" "pyproj>=3.6.0"

Confirm the GDAL backend before running a batch:

import rasterio
print(rasterio.__version__)     # e.g. 1.3.9
print(rasterio.gdal_version())  # e.g. 3.7.3

How a Fixed Grid Contract Forces Alignment

Reprojecting multiple rasters onto one shared grid contract Three differently-gridded source rasters flow through a reproject_to_grid step parameterised by a single shared contract of CRS, resolution, and bounds. All outputs emerge with the same transform, width, and height and stack into one aligned array. DEM · EPSG:4269 10 m · origin A NDVI · EPSG:32610 20 m · origin B Landcover · EPSG:5070 30 m · origin C reproject_to_grid one CRS · res · bounds contract identical transform stacks pixel-for-pixel
Independent reprojection versus one grid contract On the left three source rasters each call calculate_default_transform and end up with slightly different origins and pixel sizes, so their grids interleave. On the right the same three are warped to a single declared transform, width and height, so their cells coincide exactly and per-pixel arithmetic between them is meaningful. each raster picks its own grid 3 origins, 3 pixel sizes, no shared lattice stacking them resamples twice and blurs one contract, three warps transform width, height crs, nodata declared once, in config every cell coincides — one resample per input and the stack is reproducible next month

Because the destination geometry comes from the shared bounds and resolution — never from the individual source — every output carries the exact same affine, width, and height. That is what makes np.stack and per-pixel band math safe.

Recipe: reproject_to_grid with rasterio.warp

The function below reprojects one source into a caller-supplied grid contract. Run it once per raster with the same dst_crs, resolution, and bounds, and every output will register pixel-for-pixel.

import logging
from pathlib import Path

import rasterio
from rasterio.crs import CRS
from rasterio.transform import Affine, from_origin
from rasterio.warp import Resampling, calculate_default_transform, reproject

logger = logging.getLogger(__name__)


def reproject_to_grid(
    src_path: Path,
    out_path: Path,
    dst_crs: str,
    bounds: tuple[float, float, float, float],
    resolution: float | None = None,
    resampling: Resampling = Resampling.bilinear,
    dst_nodata: float | None = None,
    num_threads: int = 4,
) -> Affine:
    """
    Reproject one raster onto a fixed grid contract so that any set of rasters
    processed with the same (dst_crs, bounds, resolution) stacks pixel-for-pixel.

    Parameters
    ----------
    src_path:    Source raster on disk.
    out_path:    Destination GeoTIFF path.
    dst_crs:     Target CRS as an authority string, e.g. "EPSG:5070".
    bounds:      (left, bottom, right, top) of the shared grid in dst_crs units.
    resolution:  Square pixel size in dst_crs units. If None, a native-scale
                 value is derived from the source via calculate_default_transform.
    resampling:  Kernel — use Resampling.nearest for categorical rasters.
    dst_nodata:  Output nodata sentinel; falls back to the source nodata.
    num_threads: GDAL warp threads.

    Returns
    -------
    The destination Affine transform, identical for every source sharing the
    same contract.
    """
    target_crs = CRS.from_user_input(dst_crs)
    left, bottom, right, top = bounds

    with rasterio.open(src_path) as src:
        if src.crs is None:
            raise ValueError(f"{src_path}: missing CRS — normalise CRS before gridding")

        # Derive a sensible pixel size only when the caller leaves it unset.
        if resolution is None:
            _, ref_width, _ = calculate_default_transform(
                src.crs, target_crs, src.width, src.height, *src.bounds,
            )
            resolution = (right - left) / ref_width
            logger.info("Auto-selected resolution %.4f for %s", resolution, src_path)

        # Deterministic grid: width, height, and origin come from the contract
        # alone — never from this source — so all outputs share one transform.
        width = int(round((right - left) / resolution))
        height = int(round((top - bottom) / resolution))
        dst_transform = from_origin(left, top, resolution, resolution)

        out_nodata = dst_nodata if dst_nodata is not None else src.nodata
        profile = src.profile.copy()
        profile.update(
            crs=target_crs,
            transform=dst_transform,
            width=width,
            height=height,
            nodata=out_nodata,
            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:
            for band_idx in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, band_idx),
                    destination=rasterio.band(dst, band_idx),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst_transform,
                    dst_crs=target_crs,
                    resampling=resampling,
                    src_nodata=src.nodata,
                    dst_nodata=out_nodata,
                    num_threads=num_threads,
                )

    logger.info(
        "Reprojected %s -> %s (%dx%d @ %.3f)",
        src_path, out_path, width, height, resolution,
    )
    return dst_transform

Key Implementation Notes

  • The grid is built from bounds, not from the source. Computing width, height, and the affine from the contract with from_origin is the single decision that guarantees alignment. Deriving destination geometry per source — the intuitive move — is exactly what reintroduces sub-pixel drift.
The four assertions that prove a warp worked Four checks run after every reprojection. The CRS must equal the contract CRS exactly. The affine must match within a small tolerance rather than by equality, because floating-point arithmetic makes exact comparison brittle. The array shape must match the declared width and height. The valid-pixel count must not increase, which catches nodata bleeding into the data area. assert after the warp, not before — the contract is only real if it is checked 1 · dst.crs == contract.crs — exact equality, no ballpark match 2 · affine within 1e-6 — compare with a tolerance, never with == 3 · array.shape == (height, width) — an off-by-one row is a silent shift 4 · valid pixels did not grow — growth means nodata was interpolated in

  • calculate_default_transform is used only to suggest a resolution. It projects the source’s own bounds and would produce a different origin for each raster if used to build the grid. Restricting it to a resolution hint keeps its convenience without letting it dictate the origin.

  • from_origin(left, top, res, res) anchors on the top-left corner. GeoTIFF affines are north-up, so the y pixel size is negative internally; from_origin handles the sign. Passing bottom where top belongs silently flips the raster vertically.

  • Explicit src_nodata and dst_nodata prevent edge halos. Continuous kernels blend valid pixels with the nodata sentinel near masked boundaries. Declaring nodata on both sides lets GDAL exclude those pixels from interpolation instead of averaging them in.

  • Kernel choice must match the data type. Resampling.bilinear is a safe default for continuous surfaces but corrupts class codes on categorical rasters; pick the kernel deliberately per source, as covered in the companion guide on choosing resampling methods for categorical and continuous rasters.

  • Verify before you stack. After reprojecting a set, assert that every output’s transform, width, and height are identical. A cheap equality check turns a silent misregistration into a loud, early failure.

Troubleshooting Common Common-Grid Failures

Failure Root Cause Fix
ValueError on np.stack of outputs Rasters reprojected independently, differing by a row or column Build the grid from one shared bounds/resolution, not per source
Outputs share a CRS but misregister in band math Origins rounded independently by per-source calculate_default_transform Derive the transform with from_origin from the contract
Nodata halo around valid data Continuous kernel averaging valid pixels with the nodata sentinel Pass explicit src_nodata and dst_nodata to reproject
Raster written upside down bottom passed where top is expected in from_origin Anchor on the top-left corner: from_origin(left, top, res, res)
Fractional class codes in categorical output Resampling.bilinear used on an integer land-cover raster Pass resampling=Resampling.nearest for categorical inputs

Integration Note

reproject_to_grid slots into the reproject-and-resample stage described in Raster Alignment & Resampling Techniques, immediately after coordinate systems have been unified. Run CRS normalization across mixed datasets first so the datum transformation is already resolved when the warp runs, then feed each aligned output to the per-band routing described in aligning raster bands with rasterio and affine transforms when different bands need different kernels.