xarray vs Rasterio for Time-Series Rasters

Once a directory holds 120 scenes of the same place, the question stops being "how do I read a GeoTIFF" and becomes "who keeps track of which file is which date". This guide is for analysts and data engineers deciding whether a labelled cube is worth its overhead for stacked raster work; it sits under xarray & rioxarray Raster Cubes in Mastering Core Geospatial Python Libraries, and assumes you already know how to pull bands out of a single file the way reading multi-band TIFFs with Rasterio describes. The worked example throughout is one job — a monthly mean NDVI composite from 120 Sentinel-2 derived scenes — written both ways.

Why This Approach / What Goes Wrong

Rasterio's model stops at the edge of a file. A dataset object gives you bands, a transform, a CRS and windows into that raster, and it does those things with almost no abstraction between your code and GDAL. Nothing in the API spans files, so the moment your analysis crosses scenes you become the index: parsing dates out of filenames, grouping them, allocating accumulators, dividing sums by counts, and remembering that axis 0 of your stack means time this week and band last week. That bookkeeping is not hard, it is just unpaid and untested, and it lives in every script that touches the archive.

Two concrete failures follow. The first is memory. The obvious loop reads each scene fully and stacks the result — but a 10 m Sentinel-2 tile is 10980 × 10980 pixels, 482 MB as float32, so 120 of them is 57.8 GB and np.stack dies on any normal machine. The fix (accumulate over blocks instead of materializing) is correct and cheap, but it is another thirty lines of hand-written state.

The second failure is quieter and worse. A NumPy array from src.read() carries no coordinates, so if one scene is a half-pixel off, comes from a neighbouring tile, or was resampled to 20 m and reprojected, it stacks perfectly well and the composite is silently wrong. xarray cannot make that data good, but it makes the problem loud: arrays are aligned on coordinate labels, so mismatched grids either raise under join="exact" or produce visible NaN instead of plausible nonsense. That guarantee — not raw speed — is the main thing the cube sells. Rasterio is not the loser here either: it remains the better tool for single-scene I/O, for byte-level block control, and for writing anything you care about the encoding of, which is why the two-library comparison in Rasterio vs GDAL Python bindings still applies underneath the cube — rioxarray calls rasterio for every byte it moves.

Task-by-task verdict: a rasterio loop versus an xarray cube Six raster tasks compared across two columns. A monthly mean over 120 scenes takes about forty maintained lines in a rasterio loop but one resample call on a cube, so the cube wins. A per-pixel slope over time needs the whole 57.8 gigabyte stack in memory with rasterio but reduces chunk by chunk on a cube, so the cube wins. Opening one scene and reading four bands is a direct read in rasterio versus the same read plus scheduler overhead on a cube, so rasterio wins. Exact block-aligned input and output is byte-for-byte controllable in rasterio while cube chunks only approximate the blocks, so rasterio wins. Writing a tuned cloud optimized GeoTIFF exposes every GDAL knob in rasterio while the cube delegates and needs a lock, so rasterio wins. Scenes on mismatched grids stack silently wrong in rasterio but align on labels and raise or return NaN on a cube, so the cube wins. Same job, two tools — where each one earns its keep The task a rasterio loop eager NumPy windows an xarray cube labelled, lazily chunked Monthly mean, 120 scenes group by date, average parse dates, hold accumulators ~40 lines you now maintain .resample(time="MS").mean() the dates are already coordinates Per-pixel slope over time trend across the archive stack every date first 57.8 GB for one 10 m tile .polyfit("time", 1) reduces chunk by chunk One scene, four bands open, read, done src.read([2, 3, 4, 8]) no graph, no scheduler same read, plus a task graph labels you never index by Exact block-aligned I/O stream a huge scene block_windows(1) · Window() byte-for-byte control chunks approximate the blocks keep them multiples of 512 Write a tuned COG compression, overviews profile, predictor, overviews every GDAL knob exposed to_raster delegates downward needs a lock when threaded Scenes on mismatched grids a resampled or shifted tile stacks anyway — wrong values no exception, no warning aligns on labels → NaN or raise join="exact" makes it loud The cube wins across the time axis; rasterio wins at the file boundary — and most pipelines cross both.
Split the decision by task, not by library loyalty: three of these six rows are still rasterio's to win.

Prerequisites

python -m pip install "xarray>=2024.7" "rioxarray>=0.15" "rasterio>=1.3.9" \
  "dask[array]>=2024.7" "numpy>=1.26" "pandas>=2.2"

Install the whole set from one package manager. rioxarray and rasterio bind the same GDAL and PROJ builds, and a pip rasterio wheel dropped on top of a conda GDAL is the usual cause of an import that fails only on the deployment host.

Step-by-Step Implementation

1. The rasterio version: group by month, accumulate over blocks. This is the memory-safe form, not the naive one — nothing bigger than one block per scene is ever resident, and the grid check is explicit because nothing else will do it for you.

import re
from collections import defaultdict
from contextlib import ExitStack
from pathlib import Path

import numpy as np
import rasterio

scene_paths = sorted(Path("s2_ndvi").glob("ndvi_*.tif"))   # ndvi_20230612.tif, ...

by_month: dict[str, list[Path]] = defaultdict(list)
for path in scene_paths:
    stamp = re.search(r"(\d{8})", path.stem).group(1)
    by_month[stamp[:6]].append(path)                       # '202306' -> [Path, ...]

with rasterio.open(scene_paths[0]) as reference:
    ref_transform, ref_crs = reference.transform, reference.crs
    ref_shape = (reference.height, reference.width)
    blocks = [window for _, window in reference.block_windows(1)]

monthly_means: dict[str, np.ndarray] = {}
for month, paths in sorted(by_month.items()):
    composite = np.full(ref_shape, np.nan, dtype="float32")
    with ExitStack() as stack:
        sources = [stack.enter_context(rasterio.open(p)) for p in paths]
        for src in sources:                       # nothing else checks this for you
            if src.transform != ref_transform or src.crs != ref_crs:
                raise ValueError(f"{src.name} is not on the reference grid")
        for window in blocks:
            total = count = None
            for src in sources:
                tile = src.read(1, window=window, masked=True)
                valid = ~np.ma.getmaskarray(tile)
                if total is None:
                    total = np.zeros(tile.shape, dtype="float64")
                    count = np.zeros(tile.shape, dtype="uint16")
                total += np.where(valid, tile.filled(0.0), 0.0)
                count += valid
            row_off, col_off = int(window.row_off), int(window.col_off)
            rows, cols = total.shape
            composite[row_off:row_off + rows, col_off:col_off + cols] = np.divide(
                total, count, out=np.full(total.shape, np.nan), where=count > 0
            ).astype("float32")
    monthly_means[month] = composite

Every line is doing something necessary, and none of it is about NDVI. masked=True is what keeps the file's nodata value out of the sum; count is what makes the divisor per-pixel rather than per-month; ExitStack keeps ten dataset handles open across the block loop instead of reopening each file thousands of times.

2. The xarray version: let the filename become a coordinate, then reduce. xr.open_mfdataset opens all 120 files lazily through the same rasterio backend and calls preprocess on each one, which is where the date turns into a real time coordinate.

import re
from pathlib import Path

import pandas as pd
import rioxarray  # noqa: F401 — registers the .rio accessor
import xarray as xr

def stamp_time(ds: xr.Dataset) -> xr.Dataset:
    """Promote the acquisition date in the filename to a length-1 time dimension."""
    stem = Path(ds.encoding["source"]).stem
    date = pd.to_datetime(re.search(r"\d{8}", stem).group(), format="%Y%m%d")
    return ds.squeeze("band", drop=True).expand_dims(time=[date])

cube = xr.open_mfdataset(
    sorted(Path("s2_ndvi").glob("ndvi_*.tif")),
    engine="rasterio",
    preprocess=stamp_time,
    combine="nested",
    concat_dim="time",
    join="exact",                       # raise on any grid that does not match
    chunks={"x": 2048, "y": 2048},
    mask_and_scale=True,                # nodata -> NaN, as masked=True does above
)["band_data"].rename("ndvi")

monthly = cube.resample(time="MS").mean()      # skipna=True by default

Those two statements replace the whole accumulator loop, and they replace it with stronger behaviour: join="exact" fails on a mismatched grid, mask_and_scale=True handles nodata, and .mean() skips NaN per pixel, so the divisor is already per-pixel. The reduction is still lazy — monthly is a task graph over 120 files, and nothing has been read.

3. Understand where the memory actually goes. The cube is not magic; it is bounded because Dask never holds more than the chunks a task needs. With {"x": 2048, "y": 2048} and one time step per file, a chunk is 16 MB, and a monthly mean gathers the ten or so chunks covering that month — roughly 170 MB — before reducing. Peak resident memory is about that figure times the number of worker threads, plus the output. The hand-written block loop reaches the same bound by a different route, and the naive np.stack version reaches 57.8 GB and dies.

Resident memory while computing twelve monthly means from 120 scenes A log-scale chart of memory against scenes read. The eager path, which stacks every scene before reducing, climbs steadily from 0.48 gigabytes after one scene through 14 gigabytes at thirty scenes and crosses a 32 gigabyte container limit at about scene 67, where it raises MemoryError; a dashed continuation shows it would have needed 57.8 gigabytes. The xarray plus Dask path oscillates in a narrow sawtooth between roughly 0.7 and 1.2 gigabytes for the whole run as chunks are loaded and released. A hand-written rasterio windowed loop is flatter still at about 0.3 gigabytes. The note explains that the windowed loop is as memory-safe as the cube; what differs is the amount of bookkeeping code. Resident memory while reducing 120 scenes to 12 monthly means 64 GB 16 GB 4 GB 1 GB 0.25 GB container limit · 32 GB MemoryError at scene 67 np.stack would need 57.8 GB eager: read every scene, then reduce xarray + Dask chunks · 0.7–1.2 GB rasterio windowed loop · ≈0.3 GB 1 30 60 90 120 scenes read Log scale. The windowed loop is as memory-safe as the cube — what differs is how much bookkeeping you write.
Only the eager stack is disqualified on memory; between the other two the deciding factor is code you have to maintain, not RAM.

4. Handle grids that do not match, instead of asserting they do. The rasterio loop above raises on the first offending scene and stops, which is honest but unhelpful when a supplier delivered three tiles at 20 m. On the cube side, reproject_match warps each incoming scene onto a reference grid during preprocess, so the concat stays exact and the misfits are repaired rather than rejected.

import rioxarray
from rasterio.enums import Resampling

grid = rioxarray.open_rasterio(scene_paths[0], chunks=True).squeeze("band", drop=True)

def stamp_and_align(ds: xr.Dataset) -> xr.Dataset:
    ds = stamp_time(ds)
    if ds.rio.transform() != grid.rio.transform() or ds.rio.crs != grid.rio.crs:
        # bilinear for continuous NDVI; nearest for anything categorical
        ds = ds.rio.reproject_match(grid, resampling=Resampling.bilinear)
    return ds

Pass preprocess=stamp_and_align instead of stamp_time and the cube absorbs mixed resolutions without losing the join="exact" guarantee. Pick the reference grid deliberately — the finest scene if you are willing to interpolate, the coarsest if you refuse to invent detail — and keep the analysis CRS metric (a UTM zone or an equal-area projection, never Web Mercator, whose scale factor makes per-pixel area a function of latitude). The full grid-matching rules live in reprojecting raster cubes with reproject_match.

5. Write the results with rasterio. This is the hybrid that production pipelines converge on: the cube expresses the analysis, rasterio writes the output, because to_raster deliberately exposes only a slice of GDAL's creation options and needs a lock under threads.

import numpy as np
import rasterio
from rasterio.enums import Resampling

with rasterio.open(scene_paths[0]) as reference:
    profile = reference.profile.copy()

profile.update(
    driver="GTiff", count=1, dtype="float32", nodata=np.nan,
    tiled=True, blockxsize=512, blockysize=512,
    compress="DEFLATE", predictor=3,        # predictor 3 = floating point
    BIGTIFF="IF_SAFER",
)

for stamp in monthly["time"].values:
    layer = monthly.sel(time=stamp).compute()          # one month resident at a time
    name = np.datetime_as_string(stamp, unit="M").replace("-", "")
    with rasterio.open(f"ndvi_mean_{name}.tif", "w", **profile) as dst:
        dst.write(layer.values.astype("float32"), 1)
        dst.build_overviews([2, 4, 8, 16], Resampling.average)
        dst.update_tags(ns="rio_overview", resampling="average")
Decision tree for choosing a rasterio loop or an xarray cube A branching decision tree. The root asks how many scenes share one grid; just one leads straight to rasterio open, read and write. Many scenes leads to a second question, whether you reduce across time or band; if not, a rasterio loop per scene is the answer and it parallelizes trivially. If yes, the third question asks whether every scene has an identical transform and CRS; if not, reproject_match onto a reference grid first, then build the cube. If yes, the fourth question asks whether the whole stack fits in RAM; if it does, an eager cube without Dask has the least overhead, and if it does not, a chunked lazy cube with time set to minus one is the answer, whose peak memory is roughly one chunk per worker thread. A side panel gives the rule of thumb: one or two scenes means rasterio, three or more with maths across time means the cube pays for itself, mixed grids must be matched first, and tuned output is always written by rasterio. Which tool for this raster job? How many scenes share one grid? rasterio.open(...) one scene in, one scene out Do you reduce across time or band? rasterio loop per scene trivially parallel Identical transform and CRS on every scene? reproject_match(grid) then build the cube Does the whole stack fit in RAM? eager cube, no Dask least overhead chunks={"time": -1, "x": 2048} lazy cube · peak RAM ≈ chunk × threads Rule of thumb · 1–2 scenes → rasterio · 3+ scenes with maths   across dates → the cube pays · mixed grids → match first · tuned output → rasterio writes just one many no yes no yes yes no Scene count opens the tree; the extra dimension and the grid decide the rest.
Three questions separate the two tools: does an axis other than y and x vary, do the grids agree, and does the stack fit in memory.

Verification

The two paths must agree pixel for pixel on any month, and disagreement is the only reliable signal that one of them has a bookkeeping bug. Compare them where both are finite, because the cube writes NaN wherever every scene was cloud-masked and the loop writes NaN wherever count == 0 — the same pixels, by construction.

import numpy as np

xr_july = monthly.sel(time="2023-07-01").values.astype("float32")
loop_july = monthly_means["202307"]

assert xr_july.shape == loop_july.shape, "grids diverged — check join='exact'"

# NaN must appear in exactly the same places, then values must match
np.testing.assert_array_equal(np.isnan(xr_july), np.isnan(loop_july))

both = np.isfinite(xr_july)
np.testing.assert_allclose(xr_july[both], loop_july[both], rtol=1e-5, atol=1e-6)

print(f"July 2023: {both.sum():,} valid px, "
      f"max |diff| {np.nanmax(np.abs(xr_july - loop_july)):.2e}")
# July 2023: 118,452,301 valid px, max |diff| 4.77e-07

The residual difference is float32 accumulation order, not error: the loop sums in float64 and casts once, Dask sums per chunk in the array's own dtype. If you need bit-identical output, add .astype("float64") to the cube before reducing.

Edge Cases & Debugging

Frequently Asked Questions

How many scenes before the cube pays for itself? Roughly three, if the analysis actually crosses them. The break-even is not about volume but about whether you would otherwise be writing index bookkeeping — the moment your code parses a date from a filename or allocates an accumulator, the cube already costs less to maintain. For one scene in and one scene out, rasterio wins on every axis including startup time.

Is xarray slower than a plain rasterio loop? Per byte, no — rioxarray calls rasterio underneath, so the read path is identical. Per call, yes: building a Dask graph over 120 files costs a few hundred milliseconds and each task carries scheduler overhead, which is pure loss on a single small scene. On stacked work the cube usually wins anyway, because it reads only the chunks a reduction touches while a naive loop reads whole scenes.

Can I mix the two in one pipeline? That is the normal production shape. Express the analysis on the cube, then hand the materialized result to rasterio for the write, exactly as step 5 does — you get labelled, lazily evaluated maths plus full control over compression, block size, predictor and overviews. Going the other way, rioxarray.open_rasterio accepts any GDAL path including a /vsicurl/ URL or a VRT built by GDAL.

Does the cube really read only what it needs? Yes, when the chunks are aligned to the file's internal blocks. open_rasterio(..., chunks=True) adopts the file's own block layout; an explicit dict should stay a multiple of 512 for typical COGs, otherwise a single task straddles blocks and each block is fetched twice. A .sel() on time or space prunes whole chunks out of the graph before any request is made.

Related