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

Aligning Raster Bands with rasterio and Affine Transforms

Problem: Multi-source raster datasets rarely share the same affine transform, pixel size, or spatial origin. Stacking them into a common array without explicit reprojection produces silently wrong results — shifted pixels, edge NaN propagation, or coordinate system corruption that only surfaces later in band math or masking operations.

Why affine mismatches break spatial pipelines

An affine transform maps pixel coordinates (col, row) to real-world coordinates (x, y) via six coefficients (a, b, c, d, e, f). In a standard north-up raster, b and d are zero, a is pixel width, e is pixel height (negative), and c/f define the top-left corner origin. When datasets come from different sensors, processing chains, or tiling schemes, these coefficients rarely match.

Direct numpy operations on misaligned arrays produce:

• Silent spatial drift — pixels from two datasets represent different ground locations even when array indices match, invalidating any band ratio or difference calculation.
• Edge NaN propagation — differing extents leave unmasked fill values bleeding into overlap regions, corrupting aggregations and zonal statistics.
• Broken spatial joins — raster-to-vector or raster-to-raster intersections fail or return nonsensical results when pixel grids do not share an origin.
• Irreproducible outputs — floating-point origin differences accumulate across pipeline runs, making identical inputs produce subtly different outputs each time.

Version and environment compatibility

Install the production-tested stack:

# pip install rasterio>=1.3 numpy>=1.24

Verify GDAL driver availability before running in CI:

import rasterio
assert rasterio.__gdal_version__ >= "3.4.0", "GDAL 3.4+ required for reliable warping"

Affine transform anatomy — visual reference

The diagram below shows how the six affine coefficients map a pixel grid to a projected coordinate system, and what happens when two rasters share pixel dimensions but have offset origins.

Step-by-step alignment workflow

1. Extract metadata from all source rasters

Open each raster with rasterio.open in read mode and capture its CRS, bounding box (src.bounds), affine transform, and nodata value. Store these before closing file handles to avoid holding open file descriptors.

2. Detect and resolve CRS conflicts

Collect CRS objects into a set. If len(crs_set) > 1, you must reproject inputs to a common target before warping. When inputs use mixed EPSG codes, applying CRS Normalization Across Mixed Datasets as a prior step eliminates this check at alignment time.

3. Compute the union bounding box and target resolution

Take the spatial union of all src.bounds objects (min of minimums, max of maximums). Choose a target resolution: either user-specified or the finest pixel width found across src.transform.a values. Using the finest resolution avoids down-sampling any input below its native detail.

4. Build the unified affine transform

Use rasterio.transform.from_bounds(*union_bounds, width=W, height=H) where W and H are derived from the union extent divided by the target resolution. This produces a deterministic affine matrix reproducible on any run.

5. Reproject each band with rasterio.warp.reproject

Warp each band array into the target grid. Choose Resampling method by data type: nearest for categorical masks, bilinear for continuous surfaces, cubic for high-quality interpolation of fine-resolution imagery, average for aggregating coarser inputs.

6. Validate output shape, CRS, and bounds

After writing, re-open each output and assert src.shape == (H, W), src.crs == target_crs, and src.bounds == target_bounds within a floating-point tolerance.

Production-ready alignment function

import logging
import os
from pathlib import Path
from typing import Optional, Sequence

import numpy as np
import rasterio
from rasterio.crs import CRS
from rasterio.transform import from_bounds
from rasterio.warp import Resampling, reproject

log = logging.getLogger(__name__)


def align_raster_bands(
    input_paths: Sequence[str | Path],
    output_dir: str | Path,
    resampling: Resampling = Resampling.bilinear,
    target_crs: Optional[CRS] = None,
    target_resolution: Optional[float] = None,
) -> list[Path]:
    """
    Reproject a collection of raster files to a shared affine grid.

    Parameters
    ----------
    input_paths : sequence of file paths to source rasters (any GDAL-readable format)
    output_dir  : directory for aligned output files (created if absent)
    resampling  : rasterio.warp.Resampling method; use Resampling.nearest for
                  categorical data (land cover, masks), bilinear/cubic for continuous
    target_crs  : override CRS; required when inputs carry mixed projections
    target_resolution : target pixel size in CRS units; defaults to finest input resolution

    Returns
    -------
    list[Path] : paths to the aligned output files, in the same order as input_paths

    Raises
    ------
    ValueError : when inputs have mixed CRS and no target_crs is provided
    RuntimeError : when any output fails the post-write shape/bounds assertion
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    # --- 1. Collect metadata (no data is read yet) ---
    metas: list[dict] = []
    for path in input_paths:
        with rasterio.open(path) as src:
            metas.append(
                {
                    "path": Path(path),
                    "crs": src.crs,
                    "bounds": src.bounds,
                    "transform": src.transform,
                    "count": src.count,
                    "dtype": src.dtypes[0],
                    "nodata": src.nodata,
                }
            )

    # --- 2. Resolve CRS ---
    crs_set = {m["crs"] for m in metas}
    if len(crs_set) > 1 and target_crs is None:
        raise ValueError(
            f"Mixed CRS detected across inputs: {crs_set}. "
            "Provide target_crs explicitly or run CRS normalization first."
        )
    resolved_crs: CRS = target_crs or next(iter(crs_set))

    # --- 3. Union bounding box ---
    min_x = min(m["bounds"].left   for m in metas)
    min_y = min(m["bounds"].bottom for m in metas)
    max_x = max(m["bounds"].right  for m in metas)
    max_y = max(m["bounds"].top    for m in metas)
    union_bounds = (min_x, min_y, max_x, max_y)

    # --- 4. Target resolution (finest native pixel width) ---
    if target_resolution is None:
        target_resolution = min(abs(m["transform"].a) for m in metas)
    log.debug("Target resolution: %.6f CRS units/pixel", target_resolution)

    # --- 5. Target grid dimensions and affine transform ---
    width  = max(1, int(round((max_x - min_x) / target_resolution)))
    height = max(1, int(round((max_y - min_y) / target_resolution)))
    target_transform = from_bounds(*union_bounds, width=width, height=height)
    log.info("Target grid: %d x %d px, CRS=%s", width, height, resolved_crs)

    # --- 6. Reproject each raster band-by-band ---
    aligned: list[Path] = []
    for meta in metas:
        out_path = output_dir / f"aligned_{meta['path'].name}"
        nodata_val = meta["nodata"] if meta["nodata"] is not None else np.nan

        profile = {
            "driver": "GTiff",
            "crs": resolved_crs,
            "transform": target_transform,
            "width": width,
            "height": height,
            "count": meta["count"],
            "dtype": meta["dtype"],
            "nodata": nodata_val,
            "compress": "lzw",
            "tiled": True,
            "blockxsize": 256,
            "blockysize": 256,
        }

        with rasterio.open(meta["path"]) as src, rasterio.open(out_path, "w", **profile) as dst:
            for band_idx in range(1, meta["count"] + 1):
                src_array = src.read(band_idx)
                dst_array = np.full((height, width), fill_value=nodata_val, dtype=meta["dtype"])

                reproject(
                    source=src_array,
                    destination=dst_array,
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=target_transform,
                    dst_crs=resolved_crs,
                    resampling=resampling,
                    src_nodata=nodata_val,
                    dst_nodata=nodata_val,
                )
                dst.write(dst_array, band_idx)

        # --- 7. Post-write assertion ---
        with rasterio.open(out_path) as chk:
            if chk.width != width or chk.height != height:
                raise RuntimeError(
                    f"Shape mismatch in {out_path}: expected ({height}, {width}), "
                    f"got ({chk.height}, {chk.width})"
                )

        log.info("Wrote aligned raster: %s", out_path)
        aligned.append(out_path)

    return aligned

Key implementation notes

• fill_value before reproject — initialising dst_array with np.full(..., nodata_val) ensures that any pixel outside the source extent is filled with the declared nodata value rather than zero. Zero-fill is the most common source of phantom data artefacts after alignment.

• Band-sequential reads — the function reads one band at a time (src.read(band_idx)) rather than src.read() for all bands at once. This caps peak RAM usage to 1 × width × height × dtype_bytes per source, which matters when aligning dozens of Sentinel-2 scenes in a pipeline worker.

• Integer rounding for grid dimensions — max(1, int(round(...))) prevents width or height collapsing to zero on very small extents, and uses round rather than int to avoid systematic off-by-one errors from floating-point division.

• Resampling.nearest for masks — passing resampling=Resampling.nearest is mandatory for land-cover or cloud-mask rasters. Bilinear interpolation blends discrete class integers, producing class values that do not exist in the legend.

• Tiled GeoTIFF output — the profile sets tiled=True with 256×256 block size. Tiled layouts allow rasterio.windows to read spatial subsets efficiently in downstream steps without loading the full aligned raster into memory.

• Explicit src_nodata in reproject — omitting src_nodata causes GDAL to treat edge fill pixels as valid data during interpolation, producing a halo of interpolated values around the true data extent.

Integration note

This alignment function fits immediately after the CRS normalisation step in a raster ETL sequence. Wire it into the broader Raster Alignment & Resampling Techniques workflow by calling align_raster_bands before any band-math, masking, or time-series stacking operation. If your pipeline ingests multi-source satellite data, pair this with the upstream bulk satellite imagery download step so that all acquired files pass through alignment before entering the transform stage. For pipelines that also handle coordinate precision issues in vector layers, apply CRS normalization across mixed datasets to vector inputs in the same pre-processing pass, keeping both raster and vector layers in the same projected space before any spatial join.

Troubleshooting

Related

• Raster Alignment & Resampling Techniques — parent guide covering resampling kernel selection, windowed processing, and multi-temporal stacking patterns
• CRS Normalization Across Mixed Datasets — resolve mixed EPSG codes in both raster and vector inputs before alignment
• Converting Mixed EPSG Codes to a Unified CRS — step-by-step recipe for detecting and reprojecting mismatched coordinate systems
• Automated Vector & Raster Cleaning Workflows — full pipeline context: geometry repair, attribute harmonization, and spatial deduplication