Reprojecting Raster Cubes with rioxarray reproject_match

Two rasters of the same place are only comparable when they agree on four things — CRS, pixel size, grid origin, and shape — and reprojecting each one separately gets you the first without the other three. This guide is for anyone combining a satellite time series with a thematic layer, a DEM with a land-cover map, or model output with observations; it sits under xarray & rioxarray Raster Cubes in Mastering Core Geospatial Python Libraries, and assumes you can already open a cube with rioxarray.

Two independent raster grids collapsed onto one grid by reproject_match A before-and-after pair. On the left, a fine ten metre NDVI grid in EPSG 32631 with origin 412000 by 5641000 is drawn overlapping a coarser thirty metre soil grid in EPSG 4326 with origin 4.2135 by 50.9312; the two rulings cross each other at every cell, so no x or y coordinate label is shared and xarray has nothing to align on. On the right, after calling reproject_match, a single grid remains: the soil values have been resampled into the NDVI cells, and the CRS, affine transform and shape all come from the NDVI array while the x and y coordinate vectors are copied across verbatim rather than recomputed. reproject_match copies three properties: CRS, transform, shape before · two independent grids ndvi · 10 m · 32631 soil · 30 m · 4326 different origin, pixel size and CRS 412000 / 5641000 vs 4.2135 / 50.9312 no shared x/y label — nothing to align on .rio reproject _match after · one grid, one set of labels one shared grid soil values resampled into the ndvi cells crs · transform · shape from ndvi x / y vectors copied verbatim
A shared CRS is not alignment: only a shared origin, pixel size and shape give the two arrays coordinate labels xarray can match on.

Why This Approach / What Goes Wrong

.rio.reproject(dst_crs) answers the question "put this array into that CRS." To do so it has to invent an output grid, and it invents it from the source alone: rasterio's calculate_default_transform forward-transforms the source bounds, then picks a resolution that keeps roughly the source pixel count. Hand the same dst_crs to two arrays with different extents or different native resolutions and you get two different default grids. Both print EPSG:32631; neither shares a single x or y label with the other.

.rio.reproject_match(other) answers a different question: "put this array onto that array's grid." It reads the CRS, the affine transform and the shape off the match array and forwards them to reproject as explicit transform= and shape= arguments, so the output grid is not derived at all — it is copied. It then does one more thing that is easy to miss and impossible to live without: it assigns the match array's own x and y coordinate vectors onto the result. Recomputing pixel centres from an affine in float64 can land a few ULPs away from the vector the target array is carrying, and xarray's label alignment is exact — a difference of 1e-9 is a mismatch, and the resulting silent all-NaN arithmetic is described under xarray & rioxarray Raster Cubes. Copying the vectors removes that entire class of bug.

Three things reproject_match deliberately does not do. It does not modify the match array — that object is read-only input. It does not choose a resampling kernel for you: the default is Resampling.nearest for both reproject and reproject_match, which is correct for classes and quietly wrong for continuous surfaces. And it knows nothing about your nodata semantics, so an unset nodata becomes a fill value you did not choose.

Prerequisites

python -m pip install "rioxarray>=0.15" "xarray>=2024.7" "rasterio>=1.3.9" \
  "pyproj>=3.6" "numpy>=1.26" "netCDF4>=1.6"

Keep rasterio and pyproj from one package manager — both bind GDAL and PROJ, and a mixed install is the usual root cause of the transformation failures catalogued in fixing pyproj CRS transformation errors.

Step-by-Step Implementation

The worked example aligns a monthly NDVI cube — dims (time, y, x), 10 m, EPSG:32631 — with a categorical soil-class raster at 30 m in EPSG:4326.

1. Print the grid signature of both arrays before deciding anything.

Four properties decide whether a warp is needed and what it will cost: CRS, resolution, origin, shape. Print them together rather than eyeballing .rio.crs alone.

import rioxarray  # registers the .rio accessor on xarray objects
import xarray as xr

ndvi = xr.open_dataarray("ndvi_monthly_2024.nc", decode_coords="all")
soil_class = rioxarray.open_rasterio("soil_class_30m.tif").squeeze("band", drop=True)

def grid_signature(da, name):
    t = da.rio.transform()
    print(f"{name:<11} {str(da.rio.crs):<11} shape={da.shape} res={da.rio.resolution()} "
          f"origin=({t.c:.4f}, {t.f:.4f}) dtype={da.dtype} nodata={da.rio.nodata}")

grid_signature(ndvi, "ndvi")
grid_signature(soil_class, "soil_class")
# ndvi        EPSG:32631  shape=(12, 1200, 1600) res=(10.0, -10.0) origin=(412000.0000, 5641000.0000) dtype=float32 nodata=nan
# soil_class  EPSG:4326   shape=(520, 610) res=(0.00027, -0.00027) origin=(4.2135, 50.9312) dtype=uint8 nodata=None

decode_coords="all" is what makes xarray promote the NetCDF grid_mapping variable to a coordinate, so .rio.crs is populated instead of None. Before warping, confirm the two footprints actually overlap — a non-overlapping pair produces an all-nodata result with no error at all. Transformer.transform_bounds densifies the edges before projecting them, which matters because the corner points alone understate a rotated or curved footprint:

from pyproj import Transformer

tf = Transformer.from_crs(soil_class.rio.crs, ndvi.rio.crs, always_xy=True)
soil_bounds = tf.transform_bounds(*soil_class.rio.bounds())   # (left, bottom, right, top)
ndvi_bounds = ndvi.rio.bounds()

overlaps = (soil_bounds[0] < ndvi_bounds[2] and soil_bounds[2] > ndvi_bounds[0]
            and soil_bounds[1] < ndvi_bounds[3] and soil_bounds[3] > ndvi_bounds[1])
print("footprints overlap:", overlaps)   # footprints overlap: True

always_xy=True is mandatory here: EPSG:4326 is defined latitude-first, so without it the transformer returns your longitudes and latitudes swapped and the overlap test silently reports False. The full axis-order model is covered in Coordinate Systems with PyProj.

2. Choose which array is the match target.

The direction of the warp is an analysis decision, not a technical one. Warping coarse up to fine is pseudo-replication — every 10 m cell inherits its 30 m parent, so nine cells carry one independent observation and any per-pixel statistic computed on them has an inflated sample size. Warping fine down to coarse discards detail but keeps the statistics honest. The rule that survives review: match to the grid the answer is reported on. Per-pixel NDVI masked by soil type reports on the 10 m grid; a mean NDVI per soil cell reports on the 30 m grid.

from rasterio.enums import Resampling

# Direction A — the answer is per 10 m pixel: bring the class layer up to the NDVI grid
soil_on_ndvi = soil_class.rio.reproject_match(ndvi, resampling=Resampling.nearest)

# Direction B — the answer is per 30 m cell: bring the NDVI cube down to the soil grid
ndvi_on_soil = ndvi.rio.reproject_match(soil_class, resampling=Resampling.average)

If the reporting unit is a polygon rather than either grid, skip the alignment entirely and go straight to zonal statistics and raster sampling — matching grids first only adds a resampling step between you and the same number.

Coverage is the other half of the direction decision. reproject_match always returns the match array's exact shape, so if the source covers only part of the target extent the remainder comes back as fill, and if the source is larger it is cropped without warning. When the source only just reaches the target edges, pass padding=True to pad the source out to the match extent before the warp rather than losing a rim of boundary cells to the edge of the source array. Neither case raises, which is why the coverage percentage belongs in the verification block below rather than in a comment.

3. Pick the resampling kernel from the data's semantics.

Resampling is a rasterio enum, and the choice depends on two things: what the values mean, and whether the target is finer or coarser than the source.

Resampling kernel chosen by data semantics and warp direction A three-by-three matrix. Rows are data semantics: continuous surfaces such as NDVI or elevation, categorical classes such as land cover or soil, and counts or totals such as population. Columns are the warp direction: target finer than the source, target coarser than the source, and the common mistake. Continuous data takes bilinear when upsampling, average when downsampling, and is damaged by the nearest default. Categorical data takes nearest when upsampling, mode when downsampling, and is destroyed by any blending kernel because blending invents class codes that do not exist. Counts take nearest plus a division by the cell ratio when upsampling and Resampling.sum when downsampling, while plain average keeps the mean but loses the total. Which kernel: read the row for meaning, the column for direction what the values mean target finer · upsample target coarser · downsample the common mistake Continuous surface NDVI · elevation · °C Resampling.bilinear blends the 4 neighbours cubic for smoother curves Resampling.average mean of contributing pixels no aliasing Resampling.nearest the default — subsamples aliases on downsample Categorical classes land cover · soil · mask Resampling.nearest codes copied, not blended dtype stays integer Resampling.mode most frequent class wins ties break arbitrarily bilinear / cubic / average invents meaningless codes class 3 + class 7 → 5 Counts & totals population · pixel counts Resampling.nearest then divide by cell ratio keeps the grand total Resampling.sum needs GDAL 3.1 or newer conserves the total Resampling.average keeps the mean, not the sum populations quietly shrink Both reproject() and reproject_match() default to Resampling.nearest — passing nothing is still a choice.
The kernel is chosen twice over: once by what a pixel value means, once by whether the target grid is finer or coarser than the source.

Resampling.sum requires GDAL 3.1 or newer and Resampling.rms requires GDAL 3.3, so guard them if your wheels are older; Resampling.gauss is an overview-only mode and raises when passed to a warp. The same kernel table governs overview generation on export, covered in resampling and overviews when writing COGs.

4. Make nodata explicit on both sides of the warp.

A warp always produces pixels the source never covered — the target footprint is rarely a subset of the source, and rotation between CRSs guarantees corners. Those pixels need a fill value. If .rio.nodata is None, rioxarray picks a dtype-dependent default rather than failing, which is how a soil raster comes back full of 255 and an NDVI cube full of 3.4e38.

How a run of nodata pixels travels through a nearest and a bilinear warp Three horizontal strips of twelve pixels. The top strip is the source at ten metres with valid NDVI values on either side of a run of three nodata cells flagged 255. The middle strip is the result of a nearest-neighbour warp: the run of three invalid cells is unchanged and the edge between valid and invalid stays crisp. The bottom strip is the result of a bilinear warp: the invalid run has grown to five cells because every output pixel whose interpolation kernel touched an invalid neighbour is itself invalid, costing one extra cell on each side, and two on each side for cubic. A footer notes that when rio nodata is None the warp still needs a fill value, so rioxarray substitutes a dtype default of 255 for uint8 or 3.4e38 for float32 unless you set one yourself with rio write_nodata. The fill value is not optional — the only question is who chooses it source · 10 m nodata = 255 0.610.580.630.60 0.570.620.590.640.60 255255255 nearest edge stays crisp 3 in, 3 out bilinear invalid area grows 3 in, 5 out one extra cell lost on each side — roughly two for cubic, because every output pixel whose kernel touched invalid data is invalid If .rio.nodata is None the warp still needs a fill for pixels the source never covered. rioxarray substitutes a dtype default — 255 for uint8, 3.4e38 for float32. Set it yourself with .rio.write_nodata().
Nearest keeps the invalid region exactly where it was; every smoothing kernel widens it, which is why an NDVI cube loses a pixel of coastline on each warp.
import numpy as np

print(soil_class.rio.nodata)      # None → the warp has no fill value to work with

soil_class = soil_class.rio.write_nodata(255)         # the product's "unclassified" code
ndvi = ndvi.rio.write_nodata(np.nan, encoded=False)   # float cube, NaN is the sentinel

soil_on_ndvi = soil_class.rio.reproject_match(
    ndvi,
    resampling=Resampling.nearest,
    nodata=255,   # forwarded to .rio.reproject(); fills everything outside the source
)

print(soil_on_ndvi.dtype, soil_on_ndvi.rio.nodata)      # uint8 255
print(int((soil_on_ndvi == 255).sum()))                 # 18422

Two dtype rules follow from this. Never open a categorical raster with masked=True — it promotes uint8 to float64 and NaN, and a float class code is no longer a class code. And never pass nodata=np.nan to an integer array; choose a sentinel outside the legend instead, because NaN has no integer representation and the warp will raise or silently truncate.

5. Keep the extra dimensions — and know the limit.

reproject_match warps a cube with one non-spatial dimension without any special handling: the time coordinate, its attributes and its order all survive, because only the y/x axes are touched. A second non-spatial dimension is the hard limit — rioxarray maps the extra axis onto GDAL's band axis, and GDAL has only one.

import xarray as xr
from rioxarray.exceptions import TooManyDimensions

print(soil_on_ndvi.dims, ndvi.dims)   # ('y', 'x') ('time', 'y', 'x') — time survived

bands = xr.concat([red, nir], dim="band")   # dims ('time', 'band', 'y', 'x')
try:
    bands.rio.reproject_match(soil_class)
except TooManyDimensions as exc:
    print(type(exc).__name__)   # TooManyDimensions — only 2D and 3D arrays are supported

# Fix: warp one slice of the extra dim at a time, then re-stack
aligned = xr.concat(
    [bands.isel(band=i).rio.reproject_match(soil_class, resampling=Resampling.average)
     for i in range(bands.sizes["band"])],
    dim="band",
).assign_coords(band=bands["band"]).transpose("time", "band", "y", "x")

print(aligned.dims)   # ('time', 'band', 'y', 'x')

That loop is also the memory-safe pattern. rasterio.warp.reproject needs the source array in RAM, so warping a Dask-backed cube materialises it; iterating over the outer dimension caps peak memory at one slice. For cubes too large for even that, use the windowed strategies in reprojecting large datasets without memory errors.

Verification

A match is proven by four facts, not by the absence of an exception. Check them in order — each failure has a different cause.

Four assertions that prove two cubes share one grid A vertical staircase of four checks, each with a branch to the right explaining the failure. First, the two CRS objects must be equal; failing means the warp target had no written CRS. Second, the last two entries of the shapes must match; failing means reproject was called instead of reproject_match, so the output grid was derived from the source. Third, the affine transforms must be equal; failing means a half-pixel offset or a different resolution. Fourth, xr.align with join equals exact must not raise; failing means the coordinate labels drift in the final decimals, typical of a hand-built grid. When all four pass, arithmetic and masking across the two cubes are meaningful. Check in this order — each failure means something different a.rio.crs == b.rio.crs both sides carry a written CRS the source never had a CRS call .rio.write_crs() before warping a.shape[-2:] == b.shape[-2:] same height and width reproject() was used, not reproject_match() the grid was derived from the source a.rio.transform() == b.rio.transform() same origin, pixel size, rotation half-pixel offset or wrong resolution print both transforms and diff the six terms xr.align(a, b, join="exact") labels identical to the last bit float drift in a hand-built grid reproject_match copies the coord vectors all four pass → masking and arithmetic across the two cubes are meaningful reproject_match makes checks 1–3 true by construction and check 4 true by copying the coordinate vectors
Run all four: the first three compare metadata, the fourth compares the labels xarray will actually align on.
import numpy as np
import xarray as xr

target, warped = ndvi, soil_on_ndvi

assert warped.rio.crs == target.rio.crs, "CRS differs — write_crs() on the source first"
assert warped.shape[-2:] == target.shape[-2:], "grid shape differs — reproject_match not used"
assert warped.rio.transform() == target.rio.transform(), "origin or pixel size differs"

np.testing.assert_array_equal(warped["x"].values, target["x"].values)
np.testing.assert_array_equal(warped["y"].values, target["y"].values)
xr.align(target, warped, join="exact")   # raises ValueError on any label mismatch

# Coverage: how much of the target grid the source actually reached
print(f"covered by the source: {float((warped != 255).mean()) * 100:.1f}%")

# The real test: arithmetic across the pair must produce data, not an empty array
peat_ndvi = target.where(warped == 12)
print(f"transform: {tuple(round(v, 3) for v in tuple(warped.rio.transform())[:6])}")
print(f"peat cells retained: {float(peat_ndvi.notnull().mean()) * 100:.1f}%")
# covered by the source: 99.0%
# transform: (10.0, 0.0, 412000.0, 0.0, -10.0, 5641000.0)
# peat cells retained: 7.4%

Affine equality is exact float comparison, which is safe here only because reproject_match passes the target's own transform through unchanged. If you build a grid by hand with reproject(transform=..., shape=...), compare with np.allclose(tuple(a.rio.transform())[:6], tuple(b.rio.transform())[:6]) instead.

Edge Cases & Debugging

Frequently Asked Questions

When is plain .rio.reproject() the right call? When there is no second array to match. Converting a single scene to a metric CRS for measurement, or to Web Mercator at the tile-rendering boundary, is exactly what reproject is for — and resolution=, shape= and transform= let you pin the output grid when you know it in advance. The moment a second array is involved in the same expression, reproject_match is strictly safer, because it removes the chance that two independently derived grids differ.

Does reproject_match modify the array I pass as the match? No. The match array is read-only input: rioxarray pulls its CRS, transform and shape, and its x/y coordinate vectors, then returns a new array warped from the caller. Nothing is written back, and the match array is not loaded into memory beyond its coordinates. If you need both directions, call it twice with the roles reversed.

Which array should be the match target when the resolutions differ? Match to the grid your result is reported on. Upsampling a coarse layer to a fine grid replicates values and inflates the apparent sample size of any per-pixel statistic; downsampling with an aggregating kernel (average for continuous, mode for classes, sum for counts) throws away detail but keeps the numbers defensible. If the answer is reported per polygon rather than per pixel, neither direction is right — sample the raster at the polygons directly.

Can I match a whole Dataset, and does the dtype follow the match array? Dataset.rio.reproject_match works and warps every data variable onto the match grid, but a single resampling argument applies to all of them — so split a Dataset that mixes reflectance and a class mask, and warp each variable with its own kernel. Dtype and nodata come from the source variable, never from the match array, which is why a uint8 class layer stays uint8 after matching a float32 cube.

Related