Aligning Two Rasters with reproject_match

To put one raster onto another’s exact grid, use rio.reproject_match and pass a resampling method that suits the data class:

import rioxarray
from rasterio.enums import Resampling

reference = rioxarray.open_rasterio("B04_10m.tif", masked=True).squeeze(drop=True)
scl = rioxarray.open_rasterio("SCL_20m.tif").squeeze(drop=True)

scl_10m = scl.rio.reproject_match(reference, resampling=Resampling.nearest)
assert scl_10m.rio.transform() == reference.rio.transform()

This is the reliable way to satisfy the alignment requirement running through Handling Pixel Resolution and Scaling.

Why This Arises in Remote Sensing Workflows

Any operation that combines two rasters — band math, masking, differencing, stacking — requires them to describe the same ground at the same pixel positions. Satellite data rarely arrives that way. Sentinel-2 delivers 10 m, 20 m and 60 m bands; Landsat delivers 15 m, 30 m and resampled thermal; a DEM comes at whatever resolution the national mapping agency chose; a land-cover layer arrives in a different CRS entirely.

The dangerous part is that array libraries will happily combine two arrays of the same shape regardless of what they describe. If a 20 m band has been upsampled to the same shape as a 10 m band by a route that got the origin wrong, every subsequent operation is computed on mismatched ground and nothing raises. The error is geometric, and value-based checks are blind to it.

reproject_match closes that gap by construction: it takes the destination grid from an existing array rather than deriving one, so the only way to get a misaligned output is to pass the wrong reference.

Derived grid versus copied grid Reprojecting with a derived transform produces a grid whose origin depends on rounding, so it can sit half a pixel away from the reference. Copying the reference transform guarantees the output shares its origin, resolution, size and CRS, which is what pixel-for-pixel operations require. Where the destination grid comes from source: SCL 20 m EPSG:32636 reproject(crs, resolution) origin derived from bounds and rounded 10 m grid, origin 499985 reference origin is 499980 → half a pixel out, forever reproject_match(reference) copies crs, transform, width and height identical grid to the reference transforms compare equal safe to combine

Environment & Setup

Package Version Why
rioxarray ≥0.15 rio.reproject_match, rio.transform, rio.write_nodata
xarray ≥2023.6 Labelled arrays
rasterio ≥1.3.0 Resampling enum and the warp itself 
pip install "rioxarray>=0.15" "xarray>=2023.6" "rasterio>=1.3.0"

Complete Working Example

This helper aligns any number of source arrays onto one reference, choosing the resampling per array so a class layer never picks up an interpolating kernel by accident.

What reproject_match copies and what it leaves to you reproject_match copies the CRS, transform, width and height from the reference, which guarantees pixel alignment. It does not choose the resampling method, verify the source CRS is correct, or preserve a categorical dtype — those remain the caller's responsibility. Guaranteed by the call, versus still your job property handled by reproject_match your responsibility CRS copied from the reference that the source CRS is truthful transform and shape copied exactly — resampling kernel defaults to nearest pass the right one per data class dtype preserved where possible check class codes survived nodata carried if declared declare it before matching The right-hand column is where alignment bugs actually come from. 
import numpy as np
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling


def align_to(
    reference: xr.DataArray,
    sources: dict[str, xr.DataArray],
    *,
    categorical: set[str] = frozenset(),
) -> dict[str, xr.DataArray]:
    """Put every source on the reference grid, nearest for categorical layers."""
    aligned: dict[str, xr.DataArray] = {}
    for name, da in sources.items():
        method = Resampling.nearest if name in categorical else Resampling.bilinear
        out = da.rio.reproject_match(reference, resampling=method)

        # reproject_match writes the reference grid; verify rather than assume
        if out.rio.transform() != reference.rio.transform():
            raise RuntimeError(f"{name}: transform does not match the reference")
        if out.shape[-2:] != reference.shape[-2:]:
            raise RuntimeError(f"{name}: shape {out.shape[-2:]} != {reference.shape[-2:]}")

        out.attrs["resampling_applied"] = method.name
        aligned[name] = out
    return aligned


if __name__ == "__main__":
    red = rioxarray.open_rasterio("B04_10m.tif", masked=True).squeeze(drop=True)
    sources = {
        "swir": rioxarray.open_rasterio("B11_20m.tif", masked=True).squeeze(drop=True),
        "scl": rioxarray.open_rasterio("SCL_20m.tif").squeeze(drop=True),
        "dem": rioxarray.open_rasterio("dem_30m_utm.tif", masked=True).squeeze(drop=True),
    }

    aligned = align_to(red, sources, categorical={"scl"})

    stack = xr.Dataset({"red": red, **aligned})
    print(stack.dims, {k: v.attrs.get("resampling_applied") for k, v in aligned.items()})

Two details are load-bearing. The categorical set makes the nearest-neighbour rule explicit and reviewable instead of relying on whoever writes the next call site. And the transform assertion inside the helper converts a silent geometric error into an exception — cheap, because the comparison is between two six-tuples.

Variant Patterns

1. Aligning to a grid that does not exist yet

Sometimes there is no natural reference: you want a defined 10 m grid over a study area, and every input must land on it. Build the reference explicitly once, then match everything to it.

import numpy as np
import rioxarray
import xarray as xr
from rasterio.transform import from_origin

res = 10.0
left, top = 499980.0, 9900000.0     # snap the origin to a multiple of the resolution
width, height = 1200, 1200

grid = xr.DataArray(
    np.zeros((height, width), dtype="float32"),
    dims=("y", "x"),
    coords={
        "y": top - (np.arange(height) + 0.5) * res,
        "x": left + (np.arange(width) + 0.5) * res,
    },
)
grid.rio.write_crs("EPSG:32636", inplace=True)
grid.rio.write_transform(from_origin(left, top, res, res), inplace=True)

aligned = {k: v.rio.reproject_match(grid) for k, v in sources.items()}

Snapping the origin to a multiple of the resolution is what makes the grid reproducible: two people who build it independently get the same numbers, and tiles from different runs abut exactly.

2. Aligning across CRSs

reproject_match handles a CRS change and a resolution change in one warp, which is better than doing them in two steps because each warp resamples.

# dem is in a national grid; red is in UTM — one warp, not two
dem_utm10 = dem.rio.reproject_match(red, resampling=Resampling.bilinear)

Chaining reproject then reproject_match smooths the data twice and compounds the interpolation error, which is the error-budget argument made in Advanced Resampling and Upscaling Techniques.

3. Choosing the resampling per data class

Which method for which layer Reflectance and elevation upsample with bilinear or cubic and downsample with average. Class layers and boolean masks use nearest in both directions, or mode when downsampling if preserving the dominant class matters more than exact codes. reproject_match(resampling=…) layer upsampling downsampling reflectance bands bilinear / cubic average elevation / temperature bilinear average scene classification / land cover nearest nearest or mode boolean mask nearest max (conservative) `max` when downsampling a mask keeps a pixel masked if any contributing pixel was — the safe direction.

When Not to Align

Alignment is a warp, and a warp is lossy, so the best number of alignments is the smallest number that makes the analysis possible. Three situations are worth resisting.

The first is aligning for convenience rather than necessity. Two rasters only need a shared grid if they will be combined pixel by pixel. A per-scene statistic, a footprint check or an inventory needs nothing of the kind, and warping a 20 m band up to 10 m so that a summary looks tidy costs four times the memory for no gain.

The second is aligning the fine band down to the coarse one when the analysis is driven by the fine band. Downsampling a 10 m reflectance band to 20 m to match a class layer throws away exactly the detail the 10 m band was chosen for. Upsample the coarse layer with nearest instead: it adds no information, but it removes none either.

The third is repeated alignment inside a loop. If the same coarse layer is matched to the same reference for every date in a time series, the result is identical every time and the warp is pure waste. Align once, cache the aligned array or write it to disk, and reuse it — the same hoisting argument that governs mask reuse in Zonal Statistics and Vector–Raster Integration.

There is one situation where an extra alignment is worth its cost: producing an archive that many consumers will read. Warping once at write time, onto a documented grid, saves every downstream reader from doing it themselves and from disagreeing about how. That is a publishing decision rather than an analysis one, and it belongs with the other write-time choices in Writing and Validating Cloud-Optimized GeoTIFFs.

Verifying the Alignment

Equality of transforms is necessary but not sufficient, because it says nothing about whether the source was in the CRS it claimed. Add a semantic check:

import numpy as np

# 1. Structural: the grids are identical
assert scl_10m.rio.transform() == reference.rio.transform()
assert scl_10m.shape[-2:] == reference.shape[-2:]
assert scl_10m.rio.crs == reference.rio.crs

# 2. Semantic: a known feature lands where it should
row, col = 512, 700
x, y = reference.rio.transform() * (col + 0.5, row + 0.5)
print("map coordinate:", x, y, "class:", int(scl_10m.values[row, col]))

# 3. Categorical integrity: no invented classes
assert set(np.unique(scl_10m.values)) <= set(range(12)), "resampling invented class codes"

The third assertion is the one that catches a bilinear call on a class layer, and it costs a single pass over the aligned array.

Common Errors

MissingCRS when calling reproject_match

One of the arrays has no CRS attached — common after arithmetic that dropped the accessor’s state. Reattach with da.rio.write_crs("EPSG:32636", inplace=True) before matching.

The aligned array is full of NaN

The source does not overlap the reference extent, usually because it is in a different CRS than its metadata claims. Compare da.rio.bounds() with reference.rio.bounds() and consult Fixing EPSG Mismatches in rasterio.open.

The class layer gained values like 4.6

reproject_match used its default resampling, which interpolates. Pass Resampling.nearest explicitly for every categorical array.

Frequently Asked Questions

Q: What is the difference between reproject and reproject_match? reproject changes the CRS and lets the destination grid be derived; reproject_match copies the CRS, transform, width and height of another array, so the output is guaranteed to align pixel for pixel with it.

Q: Does reproject_match handle different resolutions? Yes. It resamples the source onto the reference grid whatever the resolution difference, which is why the resampling argument matters: upsampling a class layer with anything but nearest invents classes.

Q: Can I align a whole Dataset at once? Yes, but every variable is resampled with the same method. Split out any categorical variable and align it separately with nearest, then merge the results.