#Zarr and Chunked Array Storage for Raster

A satellite archive is not a picture. It is a four-dimensional array — time, band, row, column — that happens to be displayable as pictures, and the difference matters the moment somebody asks a question the picture format cannot answer. “Show me this scene” is a two-dimensional read, and Cloud-Optimized GeoTIFF serves it beautifully. “Show me twenty years of weekly NDVI at this field boundary” is a read along the time axis, and answering it from a directory of ten thousand GeoTIFFs means opening ten thousand files to extract one pixel from each.

Zarr exists for the second question. It stores an N-dimensional array as a grid of independently compressed chunks, each chunk a separate object in cloud storage, with a small JSON document describing the shape, the chunk grid, the dtype, and the codec pipeline. A reader computes which chunks its slice intersects and fetches exactly those. There is no file to open and no header to walk — the layout is the index.

This page sits in Compression, Chunking & Spatial Indexing because chunk shape is the raster analogue of row group sizing: both decide the granularity at which a cloud reader can skip data, and both punish a layout chosen without reference to the read pattern.

#Prerequisites

• Python 3.10+ with zarr>=3.0, xarray>=2024.7, numcodecs>=0.13, dask>=2024.8, and fsspec>=2024.6
• Object storage with a working credential path; the same store works on local disk, which is the right place to prototype a chunk shape
• A stated read pattern. Not a guess — an actual statement of what queries the store will serve. Every decision below is downstream of it, and no default is safe without it.
• Enough memory to hold several chunks per worker, since decompression is chunk-at-a-time and a 16 MB compressed chunk may expand to 100 MB or more

#Chunking Is the Whole Design

An N-dimensional array has one physical layout and many possible access patterns, and a chunk grid is the compromise you strike between them. Consider a daily temperature grid: 7,300 time steps × 1,800 rows × 3,600 columns of float32.

Chunk it as (1, 1800, 3600) — one full map per chunk — and a single-day continental read costs one 26 MB object. Excellent. But a twenty-year time series at one grid cell must fetch all 7,300 chunks, decompress 190 GB, and discard all but 29 KB of it.

Chunk it as (7300, 1, 1) and that time series costs one tiny object, while a single-day map costs 6.5 million requests. Also a catastrophe, in the other direction.

The workable answer is in between, and it comes from the read pattern rather than from a rule of thumb.

Two constraints bracket the choice. The lower bound is the request: an object-storage GET costs a round trip and a fraction of a cent regardless of size, so chunks below about a megabyte spend more on overhead than on payload. The upper bound is over-fetch: a chunk is the atomic unit of transfer and decompression, so a read that wants one pixel of a 200 MB chunk pays for 200 MB. Between one and sixteen megabytes compressed, both costs are tolerable. Inside that band, shape is what you tune.

#The Workflow

#1. Write the read pattern down before touching the data

Two sentences: “Readers pull full-resolution single-date windows of about 50 km, and twenty-year monthly series at individual field boundaries, roughly four to one in favour of the first.” That sentence determines every number below. Without it you are choosing a layout by taste.

#2. Size in bytes, then shape to the pattern

Fix the target compressed chunk size first — 4 MB is a sound starting point for most cloud raster. Divide by the expected compression ratio to get the uncompressed budget, divide by the item size to get the element count, then distribute those elements across the axes in proportion to how much of each axis a typical read consumes.

#3. Build the codec pipeline in the right order

A Zarr codec pipeline runs filters before the compressor. For numeric raster the order that matters is: an optional shuffle to group like bytes, then the compressor. Shuffle is nearly free and frequently worth more than two compression levels. Choosing the compressor itself is the subject of Blosc codec selection for Zarr geospatial arrays.

#4. Align writes to chunk boundaries

A write that partially covers a chunk must read the existing chunk, decompress it, patch it, recompress it, and put it back. Under concurrency, two workers doing that to the same chunk will lose one of the updates — there is no locking. Partition write jobs on chunk boundaries and the problem disappears entirely.

#5. Consolidate the metadata

Without consolidation, opening a store with 300 variables costs 300 sequential small GETs before any data moves, and on a high-latency link that is minutes. One consolidated object collapses it to a single request.

#Production Implementation

The function below derives a chunk shape from a declared read pattern and a byte budget, writes the store with a shuffle-plus-compressor pipeline, and consolidates the metadata. It refuses to write a shape whose chunks fall outside the workable band rather than producing a store that will be expensive to read.

# Requires: zarr>=3.0, xarray>=2024.7, numcodecs>=0.13, numpy>=1.26  (Python 3.10+)
from __future__ import annotations

import math
from dataclasses import dataclass

import numpy as np
import xarray as xr
import zarr
from numcodecs import Blosc

MIN_CHUNK_BYTES = 1 << 20        # 1 MiB — below this, request overhead dominates
MAX_CHUNK_BYTES = 16 << 20       # 16 MiB — above this, over-fetch dominates


@dataclass(frozen=True)
class ReadPattern:
    """How much of each axis a typical read consumes, as a fraction 0-1."""
    time: float
    y: float
    x: float

    def weights(self) -> tuple[float, float, float]:
        total = self.time + self.y + self.x
        if total <= 0:
            raise ValueError("read pattern must consume some part of some axis")
        return (self.time / total, self.y / total, self.x / total)


def chunk_shape_for(
    shape: tuple[int, int, int],
    itemsize: int,
    pattern: ReadPattern,
    *,
    target_bytes: int = 4 << 20,
    expected_ratio: float = 3.0,
) -> tuple[int, int, int]:
    """Distribute a byte budget across axes in proportion to the read pattern.

    An axis a reader traverses gets a long chunk edge; an axis it slices gets
    a short one. Returns an uncompressed chunk shape.
    """
    if target_bytes < MIN_CHUNK_BYTES or target_bytes > MAX_CHUNK_BYTES:
        raise ValueError(
            f"target {target_bytes} lies outside the workable 1-16 MiB band"
        )
    uncompressed_budget = target_bytes * expected_ratio
    elements = uncompressed_budget / itemsize
    wt, wy, wx = pattern.weights()

    # Solve k so that (k*wt)*(k*wy)*(k*wx) == elements, then clamp per axis.
    k = (elements / max(wt * wy * wx, 1e-12)) ** (1 / 3)
    raw = (k * wt, k * wy, k * wx)
    chunks = tuple(
        max(1, min(int(dim), int(math.floor(edge))))
        for dim, edge in zip(shape, raw)
    )
    return chunks  # type: ignore[return-value]


def write_zarr_store(
    dataset: xr.Dataset,
    target: str,
    pattern: ReadPattern,
    *,
    variable: str,
    clevel: int = 5,
) -> dict[str, object]:
    """Write a chunked, shuffled, consolidated Zarr store and report the layout."""
    if variable not in dataset:
        raise KeyError(f"{variable!r} is not in the dataset")
    array = dataset[variable]
    if array.ndim != 3:
        raise ValueError("expected a (time, y, x) array")

    chunks = chunk_shape_for(array.shape, array.dtype.itemsize, pattern)
    chunk_bytes = int(np.prod(chunks)) * array.dtype.itemsize
    if chunk_bytes > MAX_CHUNK_BYTES * 8:
        raise ValueError(f"uncompressed chunk of {chunk_bytes} bytes is unmanageable")

    # SHUFFLE regroups bytes of like significance before the compressor sees
    # them; on float32 raster it is worth more than two extra clevels.
    compressor = Blosc(cname="zstd", clevel=clevel, shuffle=Blosc.SHUFFLE)

    encoding = {variable: {"chunks": chunks, "compressor": compressor}}
    try:
        dataset.to_zarr(target, mode="w", encoding=encoding, consolidated=True)
    except (ValueError, OSError) as exc:
        raise RuntimeError(f"zarr write to {target} failed: {exc}") from exc

    store = zarr.open(target, mode="r")
    written = store[variable]
    return {
        "chunks": chunks,
        "uncompressed_chunk_bytes": chunk_bytes,
        "stored_bytes": int(written.nbytes_stored()),
        "logical_bytes": int(written.nbytes),
        "ratio": round(written.nbytes / max(written.nbytes_stored(), 1), 2),
        "chunk_count": int(np.prod([
            math.ceil(d / c) for d, c in zip(array.shape, chunks)
        ])),
    }

Verifying the layout is as important as choosing it. The check that matters is chunks touched per representative read:

# Requires: zarr>=3.0 — count the chunks a real slice would touch
import math


def chunks_touched(shape, chunks, selection) -> int:
    """How many chunks a slice intersects — the number that predicts cost."""
    n = 1
    for dim, chunk, sel in zip(shape, chunks, selection):
        start, stop, _ = sel.indices(dim)
        if stop <= start:
            return 0
        n *= math.floor((stop - 1) / chunk) - math.floor(start / chunk) + 1
    return n


# A 20-year monthly series at one grid cell, on a (24, 256, 256) chunk grid:
print(chunks_touched((7300, 1800, 3600), (24, 256, 256),
                     (slice(0, 7300), slice(900, 901), slice(1800, 1801))))
# → 305 chunks; at ~1.5 MB each that is ~450 MB transferred to return 29 KB.
# Widen the time edge to 512 and it drops to 15 chunks.

#Reference Matrix

Measured on a daily 1 km land-surface temperature cube, float32, 7,300 × 1,800 × 3,600, stored on S3 and read from compute in the same region.

The last row deserves emphasis because teams resist it. Storing the same array twice, chunked two ways, roughly doubles storage cost — and object storage is the cheapest thing in the stack. If both access patterns are hot, two stores beat one compromised store on every axis except a storage bill that was never the binding constraint. The same reasoning drives the tiering decisions in tiered storage for large spatial datasets.

#Failure Modes and Gotchas

Inheriting a chunk shape from the source files. Converting a directory of daily GeoTIFFs with default settings produces a store chunked one-day-per-chunk, because that is how the inputs were shaped. The conversion succeeds, the store is valid, and it is precisely the layout that makes time-series reads unusable — which was the reason for converting. Set the chunk shape explicitly in the encoding, every time.

Concurrent partial-chunk writes. Two workers writing different regions that share a chunk will each read-modify-write it, and the second put silently discards the first worker’s data. There is no error and no warning. Partition the write job on chunk boundaries; if the job cannot be partitioned that way, serialise the writes.

Forgetting to re-consolidate after appending. Consolidated metadata is a snapshot. Append a new time step without refreshing it and readers using the consolidated path will not see the new data, while readers bypassing it will — an inconsistency that is maddening to debug. Re-consolidate at the end of every write job.

Ignoring the uncompressed chunk size. A 4 MB compressed chunk of highly repetitive integer data may expand to 400 MB. Workers sized for the compressed figure will hit the memory ceiling under concurrency. Budget memory on the uncompressed size times the concurrency, not on the stored size.

Reaching for Zarr when a COG would do. A single two-dimensional orthophoto read by spatial window is exactly the case Cloud-Optimized GeoTIFF was designed for, and it carries an ecosystem of readers, tile servers, and viewers that a Zarr store does not. Use Zarr when a third axis is genuinely part of the query; do not use it because it is newer.

#Frequently Asked Questions

When is Zarr a better fit than Cloud-Optimized GeoTIFF?

Whenever the data has more than two dimensions that readers actually slice. A COG is superb for a single two-dimensional image read by spatial window, which is the common case for basemaps and orthophotos. Zarr wins when a third or fourth axis — time, spectral band, ensemble member, depth — is part of the query, because it can chunk along those axes so a time-series read touches a handful of objects instead of opening thousands of separate images.

How large should a Zarr chunk be?

Between roughly one and sixteen megabytes compressed. Below one megabyte the fixed cost of an object-storage request — TLS, round trip, per-request billing — dominates the useful payload, and a modest read turns into thousands of tiny GETs. Above about sixteen megabytes any read that wants a sliver of a chunk still pays for the whole chunk, so over-fetching swamps the saving. Within that band, shape matters far more than exact size.

What does the shuffle filter actually do before compression?

It regroups bytes so that similar bytes sit together. In a float32 array the first byte of every value is the sign and high exponent bits, which barely change across neighbouring pixels, while the last byte is noisy. Interleaved, the compressor sees an alternating pattern it cannot exploit; after a byte shuffle it sees a long run of near-identical high bytes followed by the noisy ones, and the run compresses. On typical scientific raster this is worth ten to forty per cent on its own, at almost no CPU cost.

Why does opening a Zarr store sometimes take longer than reading the data?

Because without consolidated metadata the client fetches one small JSON object per array and per group to discover the hierarchy, and a store with hundreds of variables therefore costs hundreds of sequential round trips before a single data byte moves. Writing consolidated metadata collapses that into one request. It is the single highest-value configuration change on a store with many variables and it costs one line at write time.

#Related

• Choosing Zarr Chunk Shapes for Time-Series Raster — the sizing arithmetic worked through on real cubes
• Blosc Codec Selection for Zarr Geospatial Arrays — shuffle filters, compressor choice, and decode throughput
• Zarr v3 vs COG for Multidimensional Data — the format decision, stated as a set of tests
• Cloud-Optimized GeoTIFF for Raster Workloads — the two-dimensional alternative and its tiling model
• Row Group Sizing Strategies for Parquet — the same granularity argument on the vector side

← Back to Compression, Chunking & Spatial Indexing