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

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

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.

• 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

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.

Related

• Raster Alignment & Resampling Techniques — parent guide covering the full align, reproject, validate, and write pipeline
• Choosing Resampling Methods for Categorical and Continuous Rasters — pick the right kernel for each raster before warping it onto the shared grid
• Aligning Raster Bands with rasterio and Affine Transforms — per-band kernel routing for multi-sensor stacks
• CRS Normalization Across Mixed Datasets — unify coordinate systems before reprojecting to a common grid