#Choosing Zarr Chunk Shapes for Time-Series Raster

A chunk shape is a claim about how the data will be read, and the cost of being wrong is measured in orders of magnitude rather than percentages. The same 190 GB temperature cube answers a twenty-year point series in 0.1 seconds or 21 minutes depending on a three-number decision made once at write time. This page is the arithmetic behind that decision, extending Zarr and chunked array storage for raster.

#Quick Reference

#The Arithmetic

Three constraints bound the choice, and they resolve in order.

The byte budget comes first. A compressed chunk below about a megabyte spends more on request overhead than payload; above about sixteen megabytes, any partial read over-fetches badly. Pick a target inside that band — 4 MB is a good default — and divide by the expected compression ratio to get an uncompressed byte budget, then by the item size to get an element budget.

The read pattern distributes it. If a typical read consumes 100% of the time axis and 0.05% of each spatial axis, the chunk should be long in time and short in space, in roughly that proportion. Expressed as weights that sum to one, the element budget E splits so that the product of the three edges equals E while their ratios match the weights.

Alignment adjusts it. Having derived edges of, say, (418, 137, 137), round them to something the data and the queries actually use: 365 for a year of daily steps, 256 for a spatial tile boundary. A chunk grid aligned to the query grid halves the reads for the most common query, which is worth more than a few per cent of byte-size optimality.

#Implementation

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

import math
from dataclasses import dataclass

MIN_CHUNK_BYTES = 1 << 20
MAX_CHUNK_BYTES = 16 << 20


@dataclass(frozen=True)
class AxisWeights:
    """Fraction of each axis a representative read consumes."""
    time: float
    y: float
    x: float

    def normalised(self) -> tuple[float, float, float]:
        total = self.time + self.y + self.x
        if total <= 0:
            raise ValueError("weights must sum to something positive")
        return self.time / total, self.y / total, self.x / total


def derive_chunks(
    shape: tuple[int, int, int],
    itemsize: int,
    weights: AxisWeights,
    *,
    target_compressed: int = 4 << 20,
    expected_ratio: float = 3.0,
    align: tuple[int, int, int] | None = None,
) -> tuple[int, int, int]:
    """Distribute a byte budget across axes, then snap to useful boundaries."""
    if not MIN_CHUNK_BYTES <= target_compressed <= MAX_CHUNK_BYTES:
        raise ValueError(f"{target_compressed} bytes is outside the 1-16 MiB band")

    elements = (target_compressed * expected_ratio) / itemsize
    wt, wy, wx = weights.normalised()
    k = (elements / max(wt * wy * wx, 1e-12)) ** (1 / 3)
    raw = (k * wt, k * wy, k * wx)

    chunks = []
    for dim, edge, step in zip(shape, raw, align or (1, 1, 1)):
        value = max(1, min(dim, int(math.floor(edge))))
        if step > 1:
            # Snap down to the nearest multiple so a query aligned to `step`
            # lands inside one chunk rather than straddling two.
            value = max(step, (value // step) * step)
            value = min(value, dim)
        chunks.append(value)
    return tuple(chunks)  # type: ignore[return-value]


def chunks_touched(
    shape: tuple[int, ...], chunks: tuple[int, ...], selection: tuple[slice, ...]
) -> int:
    """The number that predicts cost: how many chunks a slice intersects."""
    if not len(shape) == len(chunks) == len(selection):
        raise ValueError("shape, chunks and selection must have the same rank")
    total = 1
    for dim, chunk, sel in zip(shape, chunks, selection):
        start, stop, _ = sel.indices(dim)
        if stop <= start:
            return 0
        total *= math.floor((stop - 1) / chunk) - math.floor(start / chunk) + 1
    return total


def compare_shapes(
    shape: tuple[int, int, int],
    itemsize: int,
    candidates: dict[str, tuple[int, int, int]],
    queries: dict[str, tuple[slice, slice, slice]],
) -> None:
    """Score every candidate against every representative query."""
    header = f"{'shape':<20}" + "".join(f"{name:>22}" for name in queries)
    print(header)
    for label, chunks in candidates.items():
        mb = math.prod(chunks) * itemsize / 1e6
        row = f"{str(chunks):<20}"
        for selection in queries.values():
            row += f"{chunks_touched(shape, chunks, selection):>22,}"
        print(f"{row}   ({mb:.1f} MB uncompressed)")

#Validation

Score every candidate shape against every representative query before writing a byte. The table this produces settles arguments that would otherwise be settled by preference.

# Requires: the module above — score candidates against real queries
SHAPE = (7300, 1800, 3600)          # 20 years daily, 1 km global land
QUERIES = {
    "20y point series": (slice(0, 7300), slice(900, 901), slice(1800, 1801)),
    "1 day continent":  (slice(4000, 4001), slice(200, 1400), slice(400, 2800)),
    "1 year, 50 km box": (slice(3650, 4015), slice(900, 950), slice(1800, 1850)),
}
CANDIDATES = {
    "map-shaped":    (1, 1800, 3600),
    "series-shaped": (7300, 8, 8),
    "balanced":      (128, 512, 512),
    "aligned":       (365, 128, 128),
}
compare_shapes(SHAPE, 4, CANDIDATES, QUERIES)

Expected pattern: no shape wins every column, the aligned shape wins the year-window query decisively, and any candidate scoring above roughly 1,000 chunks on a query you actually run should be rejected outright — at 4 MB each that is gigabytes of transfer for a small answer.

#Edge Cases and Caveats

Chunk edges larger than the array. Requesting a 365-step time edge on a 200-step array silently produces a single chunk spanning everything, and the store then behaves like an unchunked file. Clamp every edge to the dimension, as the implementation does, and log when clamping happens.

Ragged final chunks. When a dimension is not a multiple of the chunk edge, the last chunk along that axis is partial — smaller, faster to read, and perfectly valid. It is only a problem if the ragged edge is tiny, which happens when the chunk edge is just under a divisor of the dimension. Nudging the edge to a clean divisor removes a class of confusing benchmark results.

Weights taken from a benchmark rather than from users. Synthetic benchmarks read uniformly; real users read a few hot regions repeatedly. If the access log is available, derive the weights from it. If it is not, state the assumption explicitly in the store’s metadata so a later reader knows what the layout was optimised for — the same discipline as recording a quantization grid.

#Frequently Asked Questions

How do I calculate how many chunks a query will read?

For each axis, take the slice’s start and stop, divide both by the chunk edge on that axis, floor them, subtract and add one — that is the number of chunks the slice spans on that axis. Multiply across axes for the total. The arithmetic takes a few lines of code and answers the only question that matters about a chunk shape, which is why it belongs in the pipeline rather than in a notebook nobody re-runs.

Should the time edge of a chunk match a calendar period?

Where readers slice by calendar period, yes, and it is worth some inefficiency to get it. If most queries ask for one year of daily data, a time edge of 365 or 366 means a year is one chunk deep instead of straddling two, halving the reads for the most common query. Aligning the chunk grid to the query grid is usually worth more than fine-tuning the byte size within the workable band.

What if two read patterns want opposite chunk shapes?

Write the array twice with different chunk grids and route each query to the store that suits it. Storage roughly doubles, which is the cheapest resource in the stack, and both patterns get their optimal layout instead of both getting a compromise that serves neither. This is only wasteful if one of the two patterns is rare, in which case give the common one the good layout and let the rare one over-fetch.

#Related

• Zarr and Chunked Array Storage for Raster — parent guide: the store layout and the consolidated-metadata step
• Blosc Codec Selection for Zarr Geospatial Arrays — the compression ratio this arithmetic assumes
• Zarr v3 vs COG for Multidimensional Data — whether to be chunking at all
• Row Group Sizing Strategies for Parquet — the same granularity argument for vector data

← Back to Zarr and Chunked Array Storage for Raster