rasterio vs xarray for Band Math

Every spectral index — NDVI, EVI, NDWI — is ultimately arithmetic on two or more bands, but the container you hold those bands in determines how much of the pipeline you have to write by hand. Raw rasterio reads give you a bare numpy array and an affine transform; you own alignment, NaN handling, and metadata propagation. Labelled band math operations with xarray give you coordinate-aware arrays that align themselves, mask themselves, and carry their CRS through every operation — at the cost of a heavier object model. This page is a decision framework for that choice, with head-to-head benchmarks and the interop pattern that lets you use both. For the broader raster foundations these two approaches share, see Core Raster Fundamentals & STAC Mapping.

The short version: reach for raw rasterio when you have one scene, a fixed grid, and you care about the last few milliseconds; reach for xarray when you have many scenes, larger-than-memory cubes, or bands that must be aligned and masked without hand-written bookkeeping.

---

Decision framework

The diagram below encodes the questions that actually change the answer. Work down it before writing any arithmetic.

The three questions that decide the outcome are: does the data fit in memory, are the inputs guaranteed to share a grid, and do you need the CRS and transform to survive the computation. If any answer pushes toward scale or safety, xarray earns its overhead. If all three are comfortable, raw numpy is the leaner tool.

---

Prerequisites

Install both stacks so you can benchmark them side by side. Pinning GDAL-linked packages avoids silent coordinate drift when the wheels update:

pip install "rasterio>=1.3" "numpy>=1.24" "xarray>=2023.1" \
            "rioxarray>=0.15" "dask[array]>=2023.1"

Conceptual prerequisites: you should understand affine geotransforms and NoData semantics from Extracting and Parsing Raster Metadata, and you should know how pixel resolution and scaling affects whether two bands even share a grid. Inputs here are surface reflectance bands with a known scale factor.

---

Head-to-head

The table is not symmetric on purpose: rasterio wins the narrow question of raw speed on one tile, while xarray wins the broad questions of correctness at scale and metadata survival. Most of the columns where xarray leads are places where the rasterio path requires code you would otherwise have to write, test, and maintain yourself.

---

When to choose raw rasterio

Choose raw rasterio and numpy when the workload is a single scene (or a windowed piece of one) that fits comfortably in memory, the bands are already on an identical grid, and you either do not need the metadata after the computation or you rebuild the profile explicitly. This path has the lowest per-tile latency because nothing constructs coordinate arrays or runs an alignment pass.

The benchmark below reads NIR and Red as plain arrays, computes NDVI with an explicit division guard, and times just the arithmetic. Note the manual float cast and the hand-written zero-denominator mask — that bookkeeping is the price of the speed.

import time
import numpy as np
import rasterio

def ndvi_rasterio(nir_path: str, red_path: str) -> tuple[np.ndarray, dict]:
    """Compute NDVI from two single-band rasters using raw numpy.
    Returns the float32 NDVI array and a rasterio profile ready for writing.
    """
    with rasterio.open(nir_path) as nir_src, rasterio.open(red_path) as red_src:
        # Shapes MUST match — raw numpy does no alignment
        if nir_src.shape != red_src.shape:
            raise ValueError(f"grid mismatch {nir_src.shape} vs {red_src.shape}")
        nir = nir_src.read(1).astype(np.float32)   # cast before arithmetic
        red = red_src.read(1).astype(np.float32)
        profile = nir_src.profile
        nodata = nir_src.nodata

    # Mask fill values by hand before the ratio
    if nodata is not None:
        nir[nir == nodata] = np.nan
        red[red == nodata] = np.nan

    denom = nir + red
    with np.errstate(divide="ignore", invalid="ignore"):
        ndvi = (nir - red) / denom
    ndvi[denom == 0] = np.nan          # explicit zero-denominator guard

    # Rebuild the profile ourselves — the array carries no metadata
    profile.update(dtype="float32", count=1, nodata=-9999.0)
    return ndvi.astype(np.float32), profile

t0 = time.perf_counter()
ndvi, profile = ndvi_rasterio("nir.tif", "red.tif")
print(f"rasterio path: {time.perf_counter() - t0:.3f}s, "
      f"valid mean={np.nanmean(ndvi):.4f}")

On a 10 980 x 10 980 Sentinel-2 tile that fits in RAM, this arithmetic runs in well under a second and allocates only the arrays you asked for. There is no hidden coordinate construction. The trade-off is visible in the code: every safety property — the shape check, the float cast, the NoData mask, the profile rebuild — is a line you own. If any is omitted, the failure is silent rather than loud.

---

When to choose xarray

Choose xarray when the data does not fit in memory, when you are processing many scenes as a stack, or when the inputs are not guaranteed to share a grid and you want alignment and masking handled for you. Backed by Dask, xarray reads and computes in chunks, so a 40 GB cube processes on a laptop; backed by rioxarray, the CRS and affine transform ride along as coordinates and survive every operation.

The benchmark below opens both bands lazily as chunked DataArrays, relies on label alignment instead of a manual shape check, and computes NDVI without ever hand-writing a mask. Nothing is read until .compute().

import time
import numpy as np
import xarray as xr
import rioxarray  # noqa: F401  registers the .rio accessor

def ndvi_xarray(nir_path: str, red_path: str, chunk: int = 2048) -> xr.DataArray:
    """Compute NDVI as a lazy, CRS-aware DataArray backed by Dask chunks."""
    nir = xr.open_dataarray(
        nir_path, engine="rasterio",
        chunks={"x": chunk, "y": chunk},
    ).squeeze("band", drop=True)
    red = xr.open_dataarray(
        red_path, engine="rasterio",
        chunks={"x": chunk, "y": chunk},
    ).squeeze("band", drop=True)

    # Convert fill values to NaN via the declared nodata; no manual np.where
    nir = nir.where(nir != nir.rio.nodata)
    red = red.where(red != red.rio.nodata)

    # Subtraction/division align on x/y labels automatically
    denom = nir + red
    ndvi = ((nir - red) / denom).where(denom != 0)
    ndvi = ndvi.rio.write_crs(nir.rio.crs)   # CRS survives the math
    return ndvi.rename("NDVI")

t0 = time.perf_counter()
ndvi = ndvi_xarray("nir.tif", "red.tif")
result = ndvi.compute()                       # trigger the Dask graph
print(f"xarray path: {time.perf_counter() - t0:.3f}s, "
      f"crs={result.rio.crs}, valid mean={float(result.mean()):.4f}")

On the same single tile this path pays a constant overhead — coordinate construction plus the alignment check — that raw numpy avoids, so it is measurably slower for one in-memory scene. Change the workload to a 24-date time-series cube and the ranking inverts: xarray streams each chunk through the same graph across every worker core, keeps peak memory near one chunk per thread, and never materialises the whole stack. For that shape of problem, the per-tile constant is noise against the throughput you gain, which is exactly why scaling raster processing with Dask is built on the xarray model.

---

Interop: read with rasterio, compute with xarray

The two are not mutually exclusive, and the strongest production pipelines use both. Validate headers with rasterio, then wrap the numpy array in a DataArray so the arithmetic gets label alignment and the metadata gets carried through. rioxarray is the bridge that attaches the transform and CRS as coordinates:

import numpy as np
import rasterio
import xarray as xr
import rioxarray  # noqa: F401

def rasterio_to_dataarray(path: str, band: int = 1) -> xr.DataArray:
    """Read one band with rasterio, hand it to xarray with CRS + transform intact."""
    with rasterio.open(path) as src:
        arr = src.read(band).astype(np.float32)
        transform, crs, nodata = src.transform, src.crs, src.nodata
        height, width = src.height, src.width

    # Build coordinate axes from the affine transform (pixel centres)
    xs = transform.c + transform.a * (np.arange(width) + 0.5)
    ys = transform.f + transform.e * (np.arange(height) + 0.5)
    da = xr.DataArray(arr, dims=("y", "x"), coords={"y": ys, "x": xs})
    da = da.rio.write_crs(crs).rio.write_nodata(nodata)
    return da

This gives you rasterio’s fast, explicit header validation and xarray’s safe arithmetic in one flow. Use it when you need the tight control of rasterio at the I/O edges but want the alignment and metadata guarantees of xarray in the middle — for example when downstream code computes NDVI directly from xarray DataArrays or stacks results into a time cube.

---

Verification & testing

Whichever path you pick, assert the same physical invariants. The two implementations should agree to floating-point tolerance on an identical, aligned tile — a divergence usually means one path skipped the NoData mask:

import numpy as np

def assert_ndvi_paths_agree(ndvi_np: np.ndarray, ndvi_xr) -> None:
    """Both computation paths must produce the same NDVI on aligned inputs."""
    xr_vals = ndvi_xr.compute().values if hasattr(ndvi_xr, "compute") else np.asarray(ndvi_xr)

    # Physical range: NDVI is bounded to [-1, 1]
    for name, arr in (("numpy", ndvi_np), ("xarray", xr_vals)):
        finite = arr[np.isfinite(arr)]
        assert finite.size > 0, f"{name}: no valid pixels"
        assert finite.min() >= -1.0 - 1e-6, f"{name}: min {finite.min():.4f} < -1"
        assert finite.max() <= 1.0 + 1e-6, f"{name}: max {finite.max():.4f} > 1"

    # Agreement where both are finite
    both = np.isfinite(ndvi_np) & np.isfinite(xr_vals)
    assert np.allclose(ndvi_np[both], xr_vals[both], atol=1e-5), \
        "paths disagree — check NoData masking on one side"
    print("both paths agree within tolerance")

Run this once on a small clipped tile in CI. If the assertion on agreement fires, the most likely cause is that the numpy path masked NoData with == while the xarray path used a declared nodata that did not match the stored fill value.

---

Pitfalls & troubleshooting

Silent grid misalignment in numpy — raw arrays with the same shape but different transforms subtract cleanly and produce a geometrically wrong result. numpy compares shapes, not coordinates. Always verify src_a.transform == src_b.transform before arithmetic, or move to xarray where alignment is on labels.

xarray “aligned to an empty region” — if two DataArrays share no overlapping coordinates, xarray’s automatic alignment intersects them to nothing and the result is all-NaN with size zero on an axis. This usually means the two rasters are in different CRSs or offset by a reprojection; reproject one with rioxarray before the operation rather than relying on alignment to rescue it.

Integer truncation on both paths — dividing two int16 bands without casting yields integer division: (NIR - Red) truncates and 0/0 gives NaN across the scene. Cast to float32 before any ratio in numpy, and open with a float dtype (or .astype) in xarray.

Lost CRS after xarray math — a plain arithmetic op on a DataArray can drop rio metadata depending on version. Re-attach with .rio.write_crs(...) after the computation, as the xarray benchmark above does, or the write step will emit a raster with no CRS.

Dask graph never runs — an xarray DataArray built from chunked reads is lazy; forgetting .compute() or a write call means no pixels are ever touched and your “fast” benchmark is really timing graph construction. Always trigger the graph when measuring.

---

FAQ

Is xarray slower than raw rasterio for band math?

For a single in-memory tile that fits in RAM, raw rasterio plus numpy is faster because it skips coordinate construction and alignment overhead — the gap is typically a few tens of milliseconds per tile. For multi-scene or larger-than-memory workloads, xarray backed by Dask wins on total throughput because it streams chunks and parallelises across cores, so raw speed on one tile is the wrong metric.

Can I mix rasterio and xarray in the same pipeline?

Yes, and it is the common production pattern. Read and validate headers with rasterio, hand the numpy array plus the affine transform and CRS to xarray by wrapping it in a DataArray, do the labelled math, then write back through rasterio or rioxarray. rioxarray is the bridge: it attaches the CRS and transform as coordinates so nothing is lost in translation.

Does xarray automatically align two bands before subtracting them?

Yes. xarray aligns on coordinate labels before any binary operation, so two DataArrays with slightly different extents are intersected to their common region and mismatched pixels become NaN. Raw numpy does not do this: it requires identical shapes and will broadcast or raise a ValueError, so grid misalignment can silently produce wrong results if you skip an explicit check.

---

Related

• Core Raster Fundamentals & STAC Mapping — the parent section covering the raster and STAC foundations both approaches build on
• Band Math Operations with xarray — the labelled-array patterns in depth, including multi-band datasets
• Calculating NDVI directly from xarray DataArrays — a worked example of the xarray path end to end
• Scaling Raster Processing with Dask — how the xarray model extends to larger-than-memory cubes and clusters
• Extracting and Parsing Raster Metadata — validate CRS, transform, and NoData before either path