Working with Zarr and Cloud-Native Datacubes

A Cloud-Optimized GeoTIFF is the right container for a scene: one date, a handful of bands, read by window. When the work is a cube — the same area across hundreds of dates, or dozens of variables on a common grid — a directory of COGs becomes awkward, because every time series read touches hundreds of files. Zarr stores the cube as one chunked, compressed array on object storage, readable lazily by xarray and Dask, and appendable as new dates arrive. This topic sits within Core Raster Fundamentals & STAC Mapping.

The one idea that governs everything else here is that a Zarr store is only as fast as its chunk shape matches the way it is read. Chunks shaped for maps make time series slow; chunks shaped for time series make maps slow. Choosing between them — or storing both — is the central design decision.

---

Prerequisites

pip install "xarray>=2024.1" "zarr>=2.17" "rioxarray>=0.15" "dask[array]>=2024.1" "fsspec>=2024.2" "s3fs>=2024.2"

You also need an aligned cube to start from — every date on one grid — which is what loading STAC items into an xarray cube with stackstac produces, and the chunking intuition from tuning Dask chunk sizes for raster cubes.

---

How a Zarr Store Is Laid Out

Two consequences follow directly from that layout. Reads parallelise naturally, because each chunk is an independent object and a selection touching forty chunks issues forty concurrent requests. And appends are cheap, because adding a date writes new chunk objects and updates the small metadata file — nothing already written is touched. Both are properties a directory of GeoTIFFs can only approximate.

---

Step-by-Step Workflow

Step 1 — Assemble an aligned cube

Zarr stores an array; it does not align one for you. Every date must share a CRS, transform and shape before it is stacked, or the cube’s x and y coordinates are a fiction.

import rioxarray
import xarray as xr


def load_aligned(paths_by_date: dict[str, str]) -> xr.DataArray:
    """Stack single-band rasters for several dates onto the first date's grid."""
    dates = sorted(paths_by_date)
    ref = rioxarray.open_rasterio(paths_by_date[dates[0]], masked=True).squeeze("band")
    frames = []
    for d in dates:
        da = rioxarray.open_rasterio(paths_by_date[d], masked=True).squeeze("band")
        if (da.rio.crs, da.rio.transform(), da.shape) != (ref.rio.crs, ref.rio.transform(), ref.shape):
            da = da.rio.reproject_match(ref)
        frames.append(da.expand_dims(time=[xr.coding.times.pd.Timestamp(d)]))
    cube = xr.concat(frames, dim="time")
    cube.name = "ndvi"
    return cube

For anything beyond a handful of dates, build the cube from a catalog search instead, which does the alignment as part of loading; the mechanics are in stackstac vs odc-stac for STAC-to-array.

Step 2 — Choose chunk shapes from the access pattern

This is the decision that determines whether the store is fast or slow, and it depends entirely on how the cube will be read.

# Shaped for maps: one date, a large spatial tile per chunk
map_chunks = {"time": 1, "y": 1024, "x": 1024}

# Shaped for time series: all dates, a small spatial tile per chunk
series_chunks = {"time": -1, "y": 128, "x": 128}

# A compromise when both matter roughly equally
balanced_chunks = {"time": 32, "y": 512, "x": 512}

The target size for any chunk is roughly 10–100 MB uncompressed. Smaller chunks mean too many objects and too many requests; larger ones mean every read fetches data it does not need. The trade-offs, and how to measure them, are the subject of choosing chunk shapes for Zarr raster cubes.

Step 3 — Write with georeferencing preserved

import numcodecs
import xarray as xr


def write_cube(cube: xr.DataArray, store: str, chunks: dict) -> None:
    ds = cube.to_dataset()
    ds = ds.rio.write_crs(cube.rio.crs)            # adds a spatial_ref variable
    ds = ds.chunk(chunks)

    encoding = {
        "ndvi": {
            "compressor": numcodecs.Blosc(cname="zstd", clevel=5,
                                          shuffle=numcodecs.Blosc.BITSHUFFLE),
            "dtype": "int16",
            "scale_factor": 0.0001,               # store NDVI as scaled int16
            "_FillValue": -32768,
        }
    }
    ds.to_zarr(store, mode="w", encoding=encoding, consolidated=True)

Three encoding choices are doing real work. Storing NDVI as scaled int16 rather than float32 halves the store with no meaningful loss of precision — the same packing described in packing float rasters into int16 with scale and offset. Blosc with ZSTD and bit-shuffling compresses smooth geophysical fields far better than plain DEFLATE. And consolidated=True gathers every metadata file into one object, so opening the store is one request rather than hundreds.

Step 4 — Append new dates

def append_date(new_frame: xr.DataArray, store: str) -> None:
    """Add one or more dates along time without touching existing chunks."""
    ds = new_frame.to_dataset(name="ndvi").rio.write_crs(new_frame.rio.crs)
    ds = ds.chunk({"time": 1, "y": 1024, "x": 1024})
    ds.drop_vars("spatial_ref").to_zarr(store, append_dim="time", consolidated=True)

Appending requires the new data to share the store’s grid exactly and to use compatible chunk shapes along the spatial dimensions. A mismatch raises rather than silently misaligning, which is one of Zarr’s more helpful behaviours. The operational details — idempotency, partial failure, time ordering — are covered in appending new time steps to a Zarr store.

Step 5 — Read lazily

import xarray as xr

cube = xr.open_zarr("s3://example-bucket/cubes/ndvi_cube.zarr", consolidated=True)

# Nothing is read yet. This selects chunks, then reads only those:
series = cube.ndvi.sel(x=780_050, y=9_950_050, method="nearest").load()
summer = cube.ndvi.sel(time=slice("2026-06-01", "2026-08-31")).median("time")

open_zarr returns a Dask-backed dataset, so every operation builds a graph and nothing is fetched until .load(), .compute() or a write forces it. A point selection touches only the chunks containing that pixel; a seasonal median touches every spatial chunk for the selected dates and nothing else.

---

Zarr and COGs Side by Side

The pattern that works best in practice is to keep both: scenes as COGs, because every viewer, tiler and desktop tool reads them and they are the natural unit of acquisition; and a Zarr cube derived from them for the analyses that are naturally cube-shaped — phenology, trend detection, per-pixel model features. The cube is a derived product, rebuildable from the COGs, which removes most of the anxiety about its format.

---

Sizing and Cost on Object Storage

A Zarr store turns one large array into many objects, and the number of objects is a real cost driver alongside the bytes. Every read is a request, every request is billed, and a store with millions of tiny chunks is slow to list, slow to copy and expensive to delete.

Compression changes the on-disk picture but not the logic. A 64 MB uncompressed chunk of scaled int16 NDVI typically stores as 10–20 MB with ZSTD and bit-shuffle, which is a comfortable object size for any object store and a comfortable read for any worker.

The other cost to plan for is rewriting. Changing the chunk shape of an existing store means reading and rewriting all of it, which for a multi-terabyte cube is a substantial job. That is the argument for getting the chunk shape right before the store grows — and, where both access patterns genuinely matter, for keeping two copies with different chunking rather than one compromise that serves neither well. Storage is usually cheaper than the compute spent working around a badly chunked cube, a trade-off examined in optimizing pipeline cost and performance.

---

Parameter Reference

---

Conventions That Keep a Cube Readable

A Zarr store is a generic container, and nothing in the format forces a cube to be self-describing. The conventions below are what make a store readable by someone other than its author a year later.

Name dimensions time, y and x, in that order, with x and y holding projected coordinates at pixel centres. Tools that infer georeferencing — rioxarray, the xarray plotting methods, most analysis libraries — look for exactly those names, and a cube with dimensions called t, row and col loses every one of those conveniences.

Write units and meaning into the variable attributes: units, long_name, the scale factor and offset if the data is packed, and a source attribute pointing at the catalog items the cube was built from. That last one is the difference between a cube that can be rebuilt and a cube that can only be trusted.

Keep one variable per physical quantity rather than stacking unrelated quantities into a band dimension. A dataset with ndvi, ndmi and cloud_fraction as separate variables is self-documenting and lets each carry its own encoding; a single array with a band dimension of three forces them to share a dtype and hides what each slice means.

---

Verification & Testing

Confirm three things after writing: the grid survived, the values round-trip, and the chunking is what you intended.

import numpy as np
import xarray as xr

back = xr.open_zarr(store, consolidated=True)

# 1. Georeferencing survived
assert back.rio.crs == cube.rio.crs, "CRS lost in the round trip"
assert np.allclose(back.x.values, cube.x.values) and np.allclose(back.y.values, cube.y.values)

# 2. Values round-trip within the packing precision
sample = dict(time=0, y=slice(0, 256), x=slice(0, 256))
diff = np.nanmax(np.abs(back.ndvi.isel(**sample).values - cube.isel(**sample).values))
assert diff <= 0.0001, f"packing error {diff} exceeds the scale factor"

# 3. Chunking is as designed
print(back.ndvi.encoding["chunks"], back.ndvi.data.chunksize)

The packing-precision check is worth keeping in any pipeline that writes scaled integers, because a wrong scale_factor — 0.001 instead of 0.0001, say — produces a cube that looks plausible and is quietly quantised ten times too coarsely.

---

Troubleshooting

ValueError: Specified zarr chunks would overlap multiple dask chunks

The Dask chunks and the Zarr chunks disagree, so parallel writers would collide on the same object. Call .chunk() with exactly the Zarr chunk shape before to_zarr.

Opening the store takes many seconds

Metadata is not consolidated, so the reader lists and fetches every metadata file. Write with consolidated=True, or run zarr.consolidate_metadata on an existing store.

The CRS is missing when the store is read back

Only the data variable was written, not the spatial_ref coordinate. Use rio.write_crs on the dataset before writing and keep the variable. Check with back.rio.crs immediately after the first write.

Appending raises a shape or chunk mismatch

The new dates are on a different grid or chunked differently in space. Reproject onto the store’s grid and rechunk to match before appending.

Time-series reads are very slow

The store is chunked for maps. Rechunk a copy for time series; the tooling and trade-offs are in reading Zarr datacubes from S3 with fsspec.

---

Frequently Asked Questions

Q: When should I use Zarr instead of COGs? When the natural unit of work is a cube rather than a scene — long time series at each pixel, multi-variable stacks, or analyses that slice along time as often as along space. COGs remain better for single-date imagery that will be tiled and viewed, because every web tool reads them.

Q: Does Zarr keep georeferencing? Only if you write it. Zarr stores arrays and attributes with no built-in notion of a CRS, so the coordinate arrays and a CRS attribute must be written alongside the data in a convention the readers understand. Using rioxarray’s spatial reference encoding is the pragmatic default.

Q: Why is reading a pixel time series from my Zarr store so slow? Because the chunks are shaped for map reads — one date and a large spatial tile each — so a single pixel’s series touches one chunk per date. Rechunk with a long time axis and small spatial extent if time series are the main access pattern.

Q: Can several writers append to one store at the same time? Only if each writes to disjoint chunks and something coordinates the metadata update. In practice, serialise the appends through one process or one scheduled task; concurrent appends to the same time dimension race on the array’s shape metadata.

---

Related

• Writing an xarray Datacube to Zarr — encodings, compressors and georeferencing in detail.
• Appending New Time Steps to a Zarr Store — incremental updates without rewriting.
• Choosing Chunk Shapes for Zarr Raster Cubes — the decision that governs performance.
• Reading Zarr Datacubes from S3 with fsspec — remote access, credentials and caching.
• Scaling Raster Processing with Dask — the execution engine behind every lazy read here.