Writing an xarray Datacube to Zarr

Attach the CRS, align the chunks, set the encoding, and write consolidated:

import numcodecs

ds = cube.to_dataset(name="ndvi").rio.write_crs(cube.rio.crs)
ds = ds.chunk({"time": 1, "y": 1024, "x": 1024})
ds.to_zarr("s3://example-bucket/cubes/ndvi.zarr", mode="w", consolidated=True,
           encoding={"ndvi": {"dtype": "int16", "scale_factor": 0.0001, "_FillValue": -32768,
                              "compressor": numcodecs.Blosc("zstd", 5, numcodecs.Blosc.BITSHUFFLE)}})

Each argument prevents a specific failure: a cube without a CRS, a write that fails on chunk overlap, a store twice the size it needs to be, and an open that takes seconds. This page belongs to working with Zarr and cloud-native datacubes in Core Raster Fundamentals & STAC Mapping.

What the Encoding Does to the Bytes

One cube, four encodings A three hundred date NDVI cube over one Sentinel-2 tile is 145 gigabytes as uncompressed float32, 62 gigabytes as float32 with default compression, 31 gigabytes as scaled int16 with default compression and 18 gigabytes as scaled int16 with ZSTD and bit-shuffling. The last is an eight-fold saving with no meaningful loss. Store size for 300 dates of NDVI over one tile float32, uncompressed 145 GB float32, default Blosc 62 GB scaled int16, default 31 GB scaled int16, ZSTD + shuffle 18 GB — eight times smaller Scaled int16 loses nothing NDVI can meaningfully express; bit-shuffle makes the compressor see the patterns.

Bit-shuffling is the setting most people have not heard of and the one that makes the biggest difference for raster data. Neighbouring pixels in a smooth field share their high-order bits; shuffling regroups the bytes so those shared bits sit next to each other, and the compressor finds long runs it would otherwise miss.

Environment & Setup

Package Version pin Used for
xarray >=2024.1 The dataset and to_zarr
zarr >=2.17 The store format
numcodecs >=0.12 Blosc compressor with ZSTD and shuffle
rioxarray >=0.15 rio.write_crs so georeferencing is stored
s3fs >=2024.2 Writing directly to object storage 
pip install "xarray>=2024.1" "zarr>=2.17" "numcodecs>=0.12" "rioxarray>=0.15" "s3fs>=2024.2"

Complete Working Example

import numcodecs
import numpy as np
import xarray as xr


def write_cube(cube: xr.DataArray, store: str, *,
               chunks: dict[str, int] | None = None,
               scale: float = 0.0001, fill: int = -32768,
               storage_options: dict | None = None) -> None:
    """Write a (time, y, x) DataArray as a georeferenced, packed, compressed Zarr store."""
    if cube.rio.crs is None:
        raise ValueError("cube has no CRS — attach one before writing")
    for dim in ("time", "y", "x"):
        if dim not in cube.dims:
            raise ValueError(f"expected dimension {dim!r}, got {cube.dims}")

    name = cube.name or "value"
    chunks = chunks or {"time": 1, "y": 1024, "x": 1024}

    ds = cube.to_dataset(name=name)
    ds = ds.rio.write_crs(cube.rio.crs)                   # spatial_ref coordinate
    ds[name].attrs.update(
        long_name=cube.attrs.get("long_name", name),
        units=cube.attrs.get("units", "1"),
        source=cube.attrs.get("source", ""),
    )
    ds = ds.chunk(chunks)                                 # Dask chunks == Zarr chunks

    lo = float(np.nanmin(cube.isel(time=0).values))
    hi = float(np.nanmax(cube.isel(time=0).values))
    if max(abs(lo), abs(hi)) / scale > 32767:
        raise ValueError(f"values {lo}..{hi} overflow int16 at scale {scale}")

    encoding = {
        name: {
            "dtype": "int16",
            "scale_factor": scale,
            "add_offset": 0.0,
            "_FillValue": fill,
            "compressor": numcodecs.Blosc(cname="zstd", clevel=5,
                                          shuffle=numcodecs.Blosc.BITSHUFFLE),
        },
        "time": {"units": "days since 2015-01-01", "calendar": "proleptic_gregorian"},
    }
    ds.to_zarr(store, mode="w", encoding=encoding, consolidated=True,
               storage_options=storage_options or {})

The overflow guard deserves its place. A scale factor that is too fine for the data range wraps values silently: NDVI of 0.9 at a scale of 0.00001 is 90,000, which does not fit in int16 and becomes a large negative number on disk. Checking one date’s range before writing catches it for the cost of reading one slice. The general rules are in avoiding integer overflow in index calculations.

Setting explicit time units matters for appends: if the first write lets xarray choose units, and a later append uses a different reference date, the concatenated axis is silently wrong.

Chunk Alignment Between Dask and Zarr

Two writers must never share a chunk When Dask chunks and Zarr chunks coincide, each Dask task writes exactly one Zarr chunk and tasks run in parallel safely. When they are misaligned, two Dask tasks each hold part of the same Zarr chunk and would overwrite each other's bytes, so xarray refuses the write rather than corrupting the store. aligned: one task, one chunk misaligned: tasks collide boundaries coincide shaded region written by two tasks Solid: Zarr chunk boundaries. Dashed: Dask chunk boundaries.

Calling .chunk() with exactly the Zarr chunk shape before writing is the simple fix, and it is why the example does so unconditionally. An alternative is to let the Dask chunks be integer multiples of the Zarr chunks — two Zarr chunks per Dask task, say — which is also safe and can reduce scheduler overhead on very large cubes.

If an upstream step produced irregular chunks — a concat of dates with different spatial chunking is the usual source — rechunking before the write is not optional. The rechunk itself can be expensive on a large cube, since it shuffles data between workers, and it is worth doing once at the end rather than repeatedly in the middle of a pipeline. The chunking trade-offs for computation, as opposed to storage, are covered in tuning Dask chunk sizes for raster cubes.

Writing Directly to Object Storage

Writing to s3:// works exactly like writing to a local path, with credentials passed through storage_options or picked up from the environment. Two details change behaviour noticeably at scale.

The first is parallelism: each Dask task writes its chunk as an independent object, so a cluster of workers writes a large cube in parallel with no coordination beyond the final metadata update. That is what makes Zarr practical for continental cubes, and it pairs naturally with the patterns in launching a Coiled cluster for STAC processing.

The second is atomicity, or rather its absence. A write that fails half way leaves a store with some chunks written and metadata describing the full shape; readers will see fill values where chunks are missing. For a store that others read, write to a new prefix and move a pointer once the write completes, rather than overwriting in place.

Verification

Three things to read back After writing, reopen the store and confirm the CRS matches, a sample of values round-trips within the packing precision, and the stored chunk shape is the one intended. Each check targets a failure that is silent at write time. Read it back before trusting it CRS present back.rio.crs == cube.rio.crs values round-trip max error within the scale chunks as intended encoding["chunks"] All three pass silently at write time when wrong; only a read-back exposes them. 
import numpy as np
import xarray as xr

back = xr.open_zarr(store, consolidated=True)
assert back.rio.crs == cube.rio.crs, "CRS did not survive"
sample = dict(time=0, y=slice(0, 256), x=slice(0, 256))
err = float(np.nanmax(np.abs(back["ndvi"].isel(**sample) - cube.isel(**sample))))
assert err <= 0.0001 + 1e-9, f"round-trip error {err}"
print("chunks:", back["ndvi"].encoding.get("chunks"))

The round-trip error should be at most one scale step. Anything larger means the scale factor, the offset or the fill value was applied inconsistently — and in particular, a fill value that collides with a real packed value turns legitimate pixels into missing data on read.

Common Errors

ValueError: Specified zarr chunks would overlap multiple dask chunks

Dask and Zarr chunk boundaries disagree. Call .chunk() with the target Zarr chunk shape immediately before writing.

Negative values appear where the data should be large

The scale factor overflowed int16. Use a coarser scale, or store as int32 or float32 for that variable.

The store opens slowly

Metadata is not consolidated. Rewrite with consolidated=True, or consolidate the existing store.

Time values are wrong after reading

Time units were not fixed at the first write. Set explicit units and calendar in the encoding.

Frequently Asked Questions

Q: Why store NDVI as int16 rather than float32? Because NDVI never needs more than four decimal places and int16 with a scale factor of 0.0001 represents it exactly to that precision at half the size. xarray applies the scale on read, so code using the cube sees floats and never has to know.

Q: Which compressor works best for raster cubes? Blosc with ZSTD and bit-shuffling is a strong default for smooth geophysical fields such as reflectance or indices. It compresses better than plain zlib and decompresses fast enough that reads remain network-bound rather than CPU-bound.

Q: What does consolidated metadata do? It gathers every array’s small metadata file into one object at the store root, so opening the store is a single request instead of one per array and coordinate. On object storage that turns a multi-second open into a fast one.

Q: Should the cube include a quality or mask variable? Yes, as its own variable with its own compact dtype. A uint8 cloud flag compresses to almost nothing and lets every consumer apply the same mask, rather than each one re-deriving it from the source scenes with slightly different rules.