xarray & rioxarray: Labelled Raster Cubes in Python

A GeoTIFF is a rectangle of numbers with a corner coordinate and a pixel size. That model is enough for a single scene, and Raster Data Handling with Rasterio covers it thoroughly. It stops being enough the moment you have forty scenes of the same place, or a climate file with pressure levels and forecast hours, or a Sentinel-2 stack where band, date and tile all vary independently — because then every array index means something, and tracking what it means becomes your problem. xarray solves that by attaching labels to every axis; rioxarray adds the geospatial half back, giving those labelled arrays a CRS, an affine transform, a nodata value, and reprojection. This part of Mastering Core Geospatial Python Libraries covers the combined model end to end: how a cube is structured, how to open one lazily, how to select by coordinate label instead of pixel index, and how to get results back out as a COG, a NetCDF file, or a Zarr store.

Architecture & Data Structures

xarray has exactly two containers. A DataArray is one N-dimensional NumPy (or Dask) array plus four pieces of metadata: dims, a tuple of axis names; coords, the actual label values along each axis; attrs, a free-form metadata dictionary; and name. A Dataset is a dict-like collection of DataArray variables that share a coordinate system — the natural fit for a NetCDF file holding temperature, precipitation and wind on one grid. Everything else in the library is built on those two.

The labels are not decoration. They drive alignment: when you subtract two DataArray objects, xarray matches them on coordinate values, not positions, and any label present in one but not the other becomes NaN. That behaviour is the source of both xarray's biggest convenience and its most confusing failure mode, and both are covered below.

rioxarray adds a .rio accessor to both containers. Importing the package is what registers it — you rarely call rioxarray by name after the import line. The accessor reads and writes a non-dimension coordinate named spatial_ref, which follows the CF convention for grid mappings, and derives everything geospatial from it plus the x/y coordinate vectors.

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

# A single scene: dims (band, y, x) with x/y in the file's own CRS units
scene = rioxarray.open_rasterio("ortho_2023_07.tif", masked=True, chunks=True)

print(scene.dims)             # ('band', 'y', 'x')
print(scene.coords["band"].values)   # [1 2 3 4] — band labels start at 1
print(scene.rio.crs)          # EPSG:32633
print(scene.rio.transform())  # Affine(10.0, 0.0, 412000.0, 0.0, -10.0, 5641000.0)
print(scene.rio.resolution()) # (10.0, -10.0)
print(scene.rio.nodata)       # None once masked=True has converted it to NaN

Three details in that snippet matter in practice. masked=True promotes the array to a float dtype and turns the file's nodata value into NaN, which is what you want before any arithmetic — an unmasked -9999 silently poisons a mean. chunks=True makes the array Dask-backed using the file's own internal block layout, so nothing is read yet. And the band coordinate is 1-based because GDAL bands are 1-based; scene.sel(band=1) and scene.isel(band=0) select the same layer, which trips people up exactly once.

There is a second entry point. xr.open_dataset(path, engine="rasterio") routes through the same rioxarray backend but returns a Dataset with a single variable called band_data, which is useful when you want a uniform Dataset interface across mixed GeoTIFF and NetCDF inputs. Passing band_as_variable=True to open_rasterio goes further and splits each band into its own named variable — handy when the bands are semantically different quantities rather than slices of one quantity.

Environment Configuration & Dependency Resolution

The stack is xarray for the labelled model, rioxarray for the geospatial accessor, rasterio (and therefore GDAL) for I/O, and Dask for out-of-core execution. rioxarray depends on rasterio and pyproj, so the C-library alignment rules from those pages apply unchanged: one package manager for the whole environment, never a pip rasterio wheel on top of a conda GDAL.

python -m pip install \
  "xarray>=2024.7" "rioxarray>=0.15" "rasterio>=1.3.9" "pyproj>=3.6" \
  "dask[array]>=2024.7" "netCDF4>=1.6" "zarr>=2.17" "bottleneck>=1.3"

netCDF4 and zarr are optional — install them only for the formats you actually read or write. bottleneck is worth having anyway: xarray dispatches rolling and NaN-aware reductions to it and gets a large speedup for free. For NetCDF specifically there are two engines, netcdf4 and h5netcdf; the latter is a pure-Python HDF5 reader that avoids the netCDF-C library entirely and is often the easier install on Windows.

import xarray as xr

# Run this in CI — prints xarray, pandas, dask, rasterio, netCDF4 and GDAL versions
xr.show_versions()

One environment variable is worth setting for cloud reads. GDAL lists the containing directory when it opens a file, which on object storage costs an extra request per open; GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR skips it. Combined with /vsicurl/ or /vsis3/ paths, open_rasterio reads remote COGs by range request exactly as described in Cloud-Native Geospatial Formats.

Vectorized Operations & Core Workflow

The canonical workflow is: open every scene lazily, concatenate them along a new time dimension, express the analysis as whole-cube operations, and trigger computation exactly once at the end. Nothing between the open and the trigger touches a pixel.

from pathlib import Path

import pandas as pd
import rioxarray
import xarray as xr

# One COG per acquisition date: ndvi_20230612.tif, ndvi_20230627.tif, ...
scene_paths = sorted(Path("s2_ndvi").glob("ndvi_*.tif"))
dates = pd.to_datetime([p.stem.split("_")[-1] for p in scene_paths], format="%Y%m%d")

layers = [
    rioxarray.open_rasterio(p, masked=True, chunks={"x": 1024, "y": 1024})
    .squeeze("band", drop=True)          # single-band scenes: drop the length-1 axis
    for p in scene_paths
]

# The new axis gets real timestamps as its coordinate, not 0..n
cube = xr.concat(layers, dim=pd.Index(dates, name="time"))
cube.name = "ndvi"
cube = cube.rio.write_crs(layers[0].rio.crs)

print(cube.dims)    # ('time', 'y', 'x')
print(cube.nbytes / 1e9, "GB if fully materialized")

xr.concat with a pd.Index as the dim argument is the idiom worth memorising: it creates the dimension and its coordinate in one step, so the resulting axis is a real DatetimeIndex that supports partial-string selection and resampling. Concatenating with dim="time" alone would produce an unlabelled axis and throw away the dates.

From there the analysis reads like NumPy with names instead of positions:

# Growing season only — a label slice, inclusive of both endpoints
summer = cube.sel(time=slice("2023-06-01", "2023-08-31"))

# Drop cloud-flagged pixels, then reduce to a monthly composite
clear = summer.where(summer > -0.2)
monthly = clear.resample(time="MS").mean()      # MS = month start

# Peak-of-season NDVI: still lazy, still nothing read
peak = monthly.max(dim="time")

# One trigger, one pass over the data
peak = peak.compute()

Because the reductions name their axis (dim="time"), the code stays correct if you later add a band dimension or reorder the axes. A positional axis=0 in the equivalent rasterio loop would not.

Chunk sizing is the one tuning knob that matters. Aim for chunks of roughly 50–200 MB and make the x/y chunk a whole multiple of the file's internal block size (512 for most COGs) so a single task never reads a partial block twice. Chunking time=1 suits per-scene work such as cloud masking; chunking time=-1 (the whole axis in one chunk) suits per-pixel time-series reductions, because each task then holds a complete series. Getting this backwards is the usual cause of a job that spends its life shuffling.

Raster Cube Processing Details

Label-based selection is the payoff for carrying coordinates around. .sel() indexes by coordinate value; .isel() indexes by integer position. Mixing them up is harmless with time, where a string cannot be an integer, and dangerous with x/y, where both are numbers.

july = cube.sel(time="2023-07")                  # partial-string datetime match
first_scene = cube.isel(time=0)                  # positional
nearest = cube.sel(time="2023-07-15", method="nearest")   # closest acquisition

# Spatial window in projected CRS units. NOTE the y slice runs high -> low,
# because north-up rasters have a DESCENDING y coordinate.
window = cube.sel(x=slice(414_000, 418_000), y=slice(5_639_000, 5_635_000))

That descending y is the single most common .sel bug on raster cubes: a slice written low-to-high silently returns a zero-length array rather than raising. Assert on window.sizes["y"] > 0 in any pipeline that builds slices from computed bounds.

Point sampling is where the labelled model produces genuinely less code than the alternative. Wrap the coordinates in DataArray objects that share a new dimension, and xarray does pointwise (not outer-product) selection:

import geopandas as gpd
import xarray as xr

sensors = gpd.read_file("sensors.gpkg").to_crs(cube.rio.crs)
site_ids = sensors["site_id"].to_numpy()

xs = xr.DataArray(sensors.geometry.x.to_numpy(), dims="site", coords={"site": site_ids})
ys = xr.DataArray(sensors.geometry.y.to_numpy(), dims="site", coords={"site": site_ids})

# Result dims: (time, site) — one full time series per sensor, in one read
series = cube.sel(x=xs, y=ys, method="nearest")
tidy = series.to_dataframe().reset_index()   # long-format DataFrame, ready for pandas

Temporal aggregation uses the same vocabulary as pandas because xarray borrows it. resample regrids the time axis onto a new frequency; groupby with a datetime component builds a climatology; rolling gives a moving window that keeps the original axis.

monthly = cube.resample(time="MS").mean()             # calendar-month composites
climatology = cube.groupby("time.month").mean("time") # mean per month-of-year
anomaly = cube.groupby("time.month") - climatology    # broadcast back against it
smoothed = cube.rolling(time=3, center=True).mean()   # 3-scene moving average

Note the pandas frequency aliases: "MS" is month-start and "YS" is year-start, both stable, whereas the old period-end alias "M" is deprecated in favour of "ME". Use the start-anchored ones and the code survives the pandas upgrade.

Reading NetCDF is the other half of the story, and the shape of the data changes: a Dataset with several variables, coordinates usually named lon/lat or longitude/latitude, and no spatial_ref at all. rioxarray needs the spatial dims named x and y (or declared explicitly) before the accessor works — the full walkthrough is in reading NetCDF climate data with xarray.

import xarray as xr
import rioxarray  # noqa: F401 — needed for .rio

ds = xr.open_dataset("era5_t2m_2023.nc", chunks={"time": 24}, engine="netcdf4")
print(list(ds.data_vars))       # ['t2m', 'tp']

t2m = (
    ds["t2m"]
    .rename({"longitude": "x", "latitude": "y"})
    .rio.set_spatial_dims(x_dim="x", y_dim="y")
    .rio.write_crs("EPSG:4326")   # ERA5 is on a plain geographic grid
)

For many files on one grid, xr.open_mfdataset(paths, combine="by_coords", parallel=True) opens them as one lazy cube without a manual concat — the right tool when the time axis is already encoded inside each file.

The head-to-head for the raster case specifically — including where a rasterio loop still wins — is worked through in xarray vs rasterio for time-series rasters.

CRS Alignment & Projection Pipeline

rioxarray stores the CRS in the spatial_ref coordinate as CF-compliant WKT, and every data variable carries a grid_mapping attribute pointing at it. That is why .rio.write_crs() returns a new object rather than mutating in place — it is adding a coordinate, not setting a field. Use inplace=True only when you deliberately want mutation.

from rasterio.enums import Resampling

# Declare a KNOWN source CRS on a cube that arrived without one (e.g. from NetCDF)
t2m = t2m.rio.write_crs("EPSG:4326")

# Pick a metric CRS from the data's own extent rather than hard-coding a zone
utm = cube.rio.estimate_utm_crs()          # e.g. EPSG:32633
cube_utm = cube.rio.reproject(utm, resampling=Resampling.bilinear)

print(cube_utm.rio.crs.to_epsg())          # 32633
print(cube_utm.rio.resolution())           # (10.0, -10.0)

Choose the resampling kernel by data semantics, exactly as in rasterio: bilinear or cubic for continuous surfaces such as reflectance or elevation, nearest for categorical rasters so class codes are never averaged into values that mean nothing. And never reproject to EPSG:3857 for measurement — its scale factor grows with latitude, so any area or distance computed there is wrong. Reproject to a UTM zone or an equal-area CRS for analysis and to Web Mercator only at the tile-render boundary. Axis order is not a concern inside rioxarray (it always works in x, y order), but the moment you build a pyproj.Transformer yourself to overlay points, pass always_xy=True.

Alignment between two cubes is a stricter requirement than a shared CRS. Because xarray aligns on coordinate labels, two arrays on the same CRS but with a half-pixel offset or a different resolution have no labels in common, and dem - landcover returns an array that is entirely NaN — no error, no warning. .rio.reproject_match() fixes this by warping one array onto the other's exact CRS, transform and shape in a single call.

import rioxarray
from rasterio.enums import Resampling

dem = rioxarray.open_rasterio("dem_10m.tif", masked=True, chunks=True).squeeze("band", drop=True)
landcover = rioxarray.open_rasterio("corine_100m.tif", chunks=True).squeeze("band", drop=True)

# Warp landcover onto the DEM's grid: same CRS, transform, width and height
landcover_on_dem = landcover.rio.reproject_match(dem, resampling=Resampling.nearest)

assert landcover_on_dem.rio.crs == dem.rio.crs
assert landcover_on_dem.shape == dem.shape
assert landcover_on_dem.rio.transform() == dem.rio.transform()

steep_forest = dem.where(landcover_on_dem == 23)   # now the arithmetic is meaningful

The grid-alignment rules, including how to choose the match target and what happens to nodata at the edges, are laid out in reprojecting raster cubes with reproject_match. For clipping a cube to a vector boundary — where the geometry CRS must be stated explicitly rather than assumed — see clipping rasters by vector geometry with rioxarray.

Production Export & Integration

Three output formats cover almost every downstream need, and the choice follows the shape of what you are exporting. A single 2-D result goes to a COG. A multi-variable cube that has to be readable by scientific tooling goes to NetCDF. A cube that will be read repeatedly, in slices, from object storage goes to Zarr.

# 1. Single 2-D layer -> Cloud Optimized GeoTIFF
peak = monthly.max(dim="time").rio.write_nodata(float("nan"), encoded=True)
peak.rio.to_raster(
    "ndvi_peak_2023.tif",
    driver="COG",          # GDAL 3.1+ writes tiling and overviews for you
    compress="DEFLATE",
    predictor=2,           # good for continuous float/int data
)

# 2. Multi-variable cube -> NetCDF with compression and int packing
cube_ds = cube.to_dataset()
cube_ds = cube_ds.rio.write_crs(cube.rio.crs)   # writes the CF spatial_ref variable
cube_ds.to_netcdf(
    "ndvi_cube_2023.nc",
    engine="netcdf4",
    encoding={"ndvi": {"zlib": True, "complevel": 4, "dtype": "int16",
                       "scale_factor": 0.0001, "_FillValue": -32768}},
)

# 3. Chunked cube -> Zarr for repeated partial reads from object storage
cube_ds.chunk({"time": 4, "y": 1024, "x": 1024}).to_zarr(
    "ndvi_cube_2023.zarr", mode="w", consolidated=True
)

Writing a large Dask-backed array to GeoTIFF deserves care: to_raster(..., windowed=True) streams block by block instead of materializing the whole array, and when several Dask threads write concurrently you need lock=threading.Lock() (or lock=True) so GDAL is not called re-entrantly on one dataset handle. Skipping the lock produces corrupt output rather than an exception.

Deployment checklist:

• Pin rioxarray and rasterio together; the accessor tracks rasterio's API and a mismatched pair fails on import, not at runtime.
• Call .rio.write_crs() immediately before every write — some xarray operations drop non-dimension coordinates, and a lost spatial_ref produces a GeoTIFF with no projection.
• Set nodata explicitly with .rio.write_nodata(value, encoded=True) so viewers and downstream masks render transparency correctly.
• Validate COG output with rio cogeo validate ndvi_peak_2023.tif before publishing; see resampling and overviews when writing COGs for the overview strategy.
• Keep attrs under control: run xr.set_options(keep_attrs=True) if provenance metadata must survive arithmetic, since xarray drops attrs by default.
• Publish Zarr with consolidated=True so a reader fetches one metadata object instead of one per array chunk.

Downstream, the same cube feeds the rest of the stack cleanly. Zonal summaries per polygon go through zonal statistics and raster sampling; a COG written here is directly consumable by a tile server, and the storage-layer trade-offs sit in Cloud-Native Geospatial Formats.

Windows / Platform Edge Cases & Debugging

Most failures here fall into two families: an HDF5 or GDAL library that does not match its Python wrapper, and a silent alignment problem that produces NaN instead of an error.

• AttributeError: 'DataArray' object has no attribute 'rio' — import rioxarray is missing. The accessor is registered by the import side effect, so linters that strip "unused" imports will break the file; add # noqa: F401.
• MissingSpatialDimensionError — the array's spatial axes are not called x/y. Rename them (.rename({"longitude": "x", "latitude": "y"})) or declare them with .rio.set_spatial_dims(x_dim="lon", y_dim="lat").
• Arithmetic returns all NaN — coordinate labels do not match. Print a.rio.transform() and b.rio.transform(); if they differ at all, run .rio.reproject_match() before the operation.
• A .sel slice returns zero elements — the y coordinate descends on north-up rasters, so the slice must run high to low. Assert on .sizes after any computed slice.
• OSError: [Errno -101] NetCDF: HDF error on a network share — HDF5 file locking. Set HDF5_USE_FILE_LOCKING=FALSE in the environment, or switch to engine="h5netcdf".
• Windows DLL load failure importing netCDF4 — a pip netCDF4 wheel layered over a conda HDF5. Rebuild the environment from a single conda-forge solve, or use h5netcdf, which has a much lighter dependency footprint.
• Corrupt or truncated GeoTIFF from a Dask write — concurrent writers without a lock. Pass lock=threading.Lock() and windowed=True to .rio.to_raster().
• Slow opens against S3 or Azure — set GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR and confirm the source is a real COG; a non-tiled GeoTIFF forces a full download per read.
• Memory blows up on a "lazy" pipeline — something triggered early. .values, plot(), to_dataframe() and even print() on a small array all compute; check for an accidental .compute() inside a loop.

Frequently Asked Questions

When is a cube actually better than a loop over rasterio? When any axis other than y and x varies and you need to operate across it. Compositing forty dates, computing a per-pixel trend, taking a monthly climatology, or masking one stack by another are all one-liners on a cube and forty-iteration loops with manual bookkeeping in rasterio. For a single scene in and a single scene out — reproject one file, write one COG — rasterio is leaner and has fewer moving parts.

Does open_rasterio read the whole file into memory? Only if you omit chunks. Without it the array is NumPy-backed and eager, so the data lands in RAM as soon as any value is touched. With chunks=True (or an explicit dict) the array is Dask-backed and the open reads only header metadata — shape, dtype, transform, CRS. Computation happens at .compute(), .load(), or when you write output.

open_rasterio or open_dataset — which should I use? Use rioxarray.open_rasterio for GDAL-readable rasters when you want a DataArray with band/y/x and the CRS already attached. Use xr.open_dataset for NetCDF, HDF5, GRIB and Zarr, where the file holds several named variables on a shared grid; add engine="rasterio" if you want a GeoTIFF presented as a Dataset for interface consistency. Both paths end up with the same .rio accessor.

Why did my CRS disappear after a groupby? Some xarray operations drop non-dimension coordinates, and spatial_ref is one. The fix is to re-attach it — result = result.rio.write_crs(source.rio.crs) — immediately before writing. Treat write_crs as cheap and call it defensively at every stage boundary rather than assuming it survived.

How big should chunks be? Target 50–200 MB per chunk and align the x/y chunk size to a multiple of the file's internal block size, usually 512 for COGs. Chunk time=1 for per-scene work such as cloud masking, and time=-1 for per-pixel time-series reductions so each task owns a whole series. Too-small chunks drown the scheduler in task overhead; too-large chunks exhaust memory when several run in parallel.

Can I write a whole cube to one GeoTIFF? Only if it reduces to at most three dimensions with the extra axis mapped onto bands — cube.rio.to_raster() will write a (time, y, x) array as a multi-band file, but the timestamps become band numbers and the labels are lost. If the time axis matters downstream, write NetCDF or Zarr; use GeoTIFF for the 2-D results you extracted from the cube.