Zarr vs NetCDF for Cloud-Native Coastal Pipelines

Operational intent: This page belongs to the Marine Spatial Data Fundamentals & Architecture domain and settles a single recurring decision — whether a multidimensional coastal or marine array should live on object storage as a Zarr store or as a NetCDF4/HDF5 file. That decision is not aesthetic. It fixes how bytes reach a worker, whether hundreds of Dask tasks can write the same array concurrently, and how much latency the first open_dataset call pays before any science runs. Coastal pipelines routinely stage decade-long ROMS hindcasts, hourly sea-surface-temperature stacks, ensemble surge forecasts, and gridded tide-gauge reanalyses — all four-dimensional (time, depth, lat, lon) arrays that are far too large to hold in RAM and are increasingly read directly from S3 or GCS. Choosing the wrong container here produces whole-file GETs on every read, metadata storms that stall job startup, and chunk-size mismatches that amplify I/O by an order of magnitude. This guide draws the boundary precisely, and deliberately stays in the array world: the sibling guide on NetCDF vs GeoTIFF for marine data covers the raster question of when a 2D grid belongs in a GeoTIFF instead, which is a different fork entirely.

---

Reference Comparison Table #

Both formats describe the same logical object — a labelled N-dimensional array with coordinates, units, and CF metadata. They differ in how that object is physically laid out and therefore in how it behaves on object storage. Zarr stores each compressed chunk as an independent key (an object on S3/GCS, a file on disk); NetCDF4 writes everything into one HDF5 container with an internal B-tree index. Read the table below as a decision reference, not a scoreboard: several rows favour Zarr for cloud pipelines, but the NetCDF/HDF5 ecosystem still wins on archival portability and server interoperability.

Two rows deserve emphasis for coastal work. Parallel region writes are what let a Prefect flow fan a 40-year hindcast across a Dask cluster and have every worker write its own year into one logical store without a lock — a pattern that is awkward-to-impossible with a single NetCDF file on S3. ERDDAP/THREDDS support is the counterweight: if your deliverable is an OPeNDAP endpoint that agency partners already point their tools at, NetCDF remains the path of least resistance, and you convert to Zarr only for the internal compute tier.

Memory-Constrained Python Implementation #

The audience runs this code. Both snippets below are typed, log through the logging module rather than print, and raise on invalid state instead of falling back silently. Neither loads the full array into RAM: xarray keeps the data lazy behind Dask until to_zarr/to_netcdf streams it chunk by chunk.

(a) Writing chunked Zarr to S3/GCS via fsspec #

This routine writes a Dask-backed dataset to an object store through an fsspec mapper. It fixes the chunk grid explicitly, attaches a Zarr compression codec, and writes consolidated metadata so the first downstream open pays a single GET rather than one-per-array.

"""
write_zarr.py — stream a Dask-backed coastal array to object storage as Zarr.

Dependencies:
    pip install "xarray>=2024.3" "zarr>=2.17" "fsspec>=2024.3" s3fs gcsfs dask numcodecs
"""
from __future__ import annotations

import logging
from typing import Mapping

import fsspec
import numcodecs
import xarray as xr

logger = logging.getLogger(__name__)

# Target ~5-20 MB compressed per chunk: small enough for one HTTP GET,
# large enough that per-object overhead stays negligible on S3/GCS.
DEFAULT_CHUNKS: dict[str, int] = {"time": 24, "depth": 1, "lat": 512, "lon": 512}


def write_coastal_zarr(
    ds: xr.Dataset,
    store_uri: str,
    chunks: Mapping[str, int] = DEFAULT_CHUNKS,
    storage_options: Mapping[str, object] | None = None,
) -> str:
    """
    Write a validated coastal Dataset to a Zarr store on object storage.

    Parameters
    ----------
    ds
        A CRS-assigned, CF-attributed xarray Dataset. Must already carry
        `_FillValue`/`units` on data variables and a `grid_mapping`.
    store_uri
        fsspec URI, e.g. ``s3://coastal-grids/roms_2024.zarr`` or
        ``gcs://coastal-grids/roms_2024.zarr``.
    chunks
        Analysis-aligned chunk grid. The SAME grid must be used by readers.
    storage_options
        Credentials/config forwarded to the fsspec backend.

    Returns
    -------
    The store URI on success.

    Raises
    ------
    ValueError
        If the Dataset has no data variables or a chunk dim is absent.
    """
    if not ds.data_vars:
        raise ValueError("Refusing to write an empty Zarr store: no data variables.")

    missing = set(chunks) - set(ds.dims)
    if missing:
        raise ValueError(f"Chunk dims absent from dataset: {sorted(missing)}")

    # Re-chunk lazily; this does not materialise data in RAM.
    ds = ds.chunk({dim: chunks[dim] for dim in chunks})

    # Explicit per-variable encoding. Blosc+zstd with a byte shuffle is a
    # strong default for float32 ocean fields; drop `_FillValue` collisions
    # by declaring it here rather than relying on inherited attrs.
    compressor = numcodecs.Blosc(cname="zstd", clevel=5, shuffle=numcodecs.Blosc.SHUFFLE)
    encoding: dict[str, dict[str, object]] = {
        var: {
            "compressor": compressor,
            "chunks": tuple(chunks.get(d, ds.sizes[d]) for d in ds[var].dims),
        }
        for var in ds.data_vars
    }

    mapper = fsspec.get_mapper(store_uri, **(storage_options or {}))
    logger.info("Writing Zarr store %s with chunks %s", store_uri, dict(chunks))

    ds.to_zarr(
        mapper,
        mode="w",
        encoding=encoding,
        consolidated=True,   # single .zmetadata read on first open
        compute=True,        # stream chunks now; swap to False for deferred Dask graphs
    )

    logger.info("Zarr store complete and metadata consolidated: %s", store_uri)
    return store_uri


def append_time_region(
    ds_new: xr.Dataset,
    store_uri: str,
    storage_options: Mapping[str, object] | None = None,
) -> None:
    """
    Append new time steps to an existing store without rewriting it.

    Uses ``append_dim='time'`` so concurrent forecast cycles can extend one
    logical store. Region writes (``region={'time': slice(...)}``) are the
    concurrent-safe alternative when the time axis is pre-allocated.

    Raises
    ------
    KeyError
        If the incoming dataset lacks a 'time' dimension.
    """
    if "time" not in ds_new.dims:
        raise KeyError("Append target requires a 'time' dimension on ds_new.")

    mapper = fsspec.get_mapper(store_uri, **(storage_options or {}))
    logger.info("Appending %d time steps to %s", ds_new.sizes["time"], store_uri)
    ds_new.to_zarr(mapper, mode="a", append_dim="time", consolidated=True)

(b) Reading and writing NetCDF-CF with xarray #

When the deliverable is a CF-compliant NetCDF4 file — for an ERDDAP dataset, a THREDDS catalog, or an agency archive — align the HDF5 chunk sizes to the same analysis grid you would have used for Zarr, and set per-variable compression explicitly. Reading uses Dask-backed chunks so a multi-gigabyte file never lands in RAM whole.

"""
netcdf_cf.py — CF-compliant NetCDF4 read/write for coastal arrays.

Dependencies:
    pip install "xarray>=2024.3" "netCDF4>=1.7" dask
"""
from __future__ import annotations

import logging
from pathlib import Path

import xarray as xr

logger = logging.getLogger(__name__)

HDF5_CHUNKS: dict[str, int] = {"time": 24, "depth": 1, "lat": 512, "lon": 512}


def read_netcdf_cf(path: str | Path) -> xr.Dataset:
    """
    Open a NetCDF4 file lazily with CF decoding and Dask chunks.

    Raises
    ------
    FileNotFoundError
        If the path does not exist.
    ValueError
        If no recognised geographic coordinate variables are present.
    """
    path = Path(path)
    if not path.exists():
        raise FileNotFoundError(f"NetCDF file not found: {path}")

    logger.info("Opening NetCDF4 (lazy): %s", path)
    ds = xr.open_dataset(
        path,
        engine="netcdf4",
        chunks={"time": 24, "lat": 512, "lon": 512},
        decode_cf=True,
        mask_and_scale=True,
    )

    coords = {c for c in ds.coords} & {"lat", "lon", "latitude", "longitude"}
    if not coords:
        raise ValueError(
            f"{path} has no recognised geographic coordinates; "
            "cannot treat as a CF grid."
        )
    return ds


def write_netcdf_cf(ds: xr.Dataset, path: str | Path) -> Path:
    """
    Write a CF-compliant NetCDF4 file with per-variable zlib compression and
    HDF5 chunk sizes aligned to the analysis grid.

    Raises
    ------
    ValueError
        If a data variable is missing a CF grid_mapping reference.
    """
    path = Path(path)

    for var in ds.data_vars:
        if "grid_mapping" not in ds[var].attrs:
            raise ValueError(
                f"Variable '{var}' lacks a CF 'grid_mapping' attribute; "
                "assign one before writing (see the CRS alignment workflow)."
            )

    encoding: dict[str, dict[str, object]] = {
        var: {
            "zlib": True,
            "complevel": 4,
            "shuffle": True,
            "chunksizes": tuple(
                min(HDF5_CHUNKS.get(d, ds.sizes[d]), ds.sizes[d]) for d in ds[var].dims
            ),
        }
        for var in ds.data_vars
    }

    logger.info("Writing CF NetCDF4: %s", path)
    ds.to_netcdf(path, engine="netcdf4", format="NETCDF4", encoding=encoding)
    logger.info("NetCDF4 written (%.1f MB): %s", path.stat().st_size / 1e6, path)
    return path

The symmetry is deliberate: the HDF5_CHUNKS grid matches DEFAULT_CHUNKS from the Zarr writer. Keeping one authoritative chunk grid across both formats means a dataset can be re-materialised in either container without the read-amplification penalty described below. Both writers assume the dataset already carries a correct CRS and grid_mapping; that contract is produced upstream by the CRS alignment workflow.

Validation Gates and Quality Control #

No store is published until it passes three gates. Each raises rather than warning, so a defective store cannot escape into the compute tier.

Gate 1 — CF Compliance #

Assert that every data variable carries units and a grid_mapping, and that coordinate axes are present. Absent CF attributes cause silent datum and unit assumptions in every downstream consumer.

import logging

import xarray as xr

logger = logging.getLogger(__name__)


def assert_cf_ready(ds: xr.Dataset) -> None:
    """Raise if the dataset is not minimally CF-compliant for publication."""
    for var in ds.data_vars:
        attrs = ds[var].attrs
        if "grid_mapping" not in attrs:
            raise ValueError(f"'{var}' missing grid_mapping (CRS reference).")
        if "units" not in attrs:
            raise ValueError(f"'{var}' missing 'units' attribute.")
    if not ({"lat", "lon", "latitude", "longitude"} & set(ds.coords)):
        raise ValueError("No geographic coordinate variables found.")
    logger.info("CF compliance gate passed for %d variables.", len(ds.data_vars))

Gate 2 — Chunk Alignment #

The single most common cloud performance defect is a mismatch between the chunk grid the writer used and the chunk grid the reader requests. If a reader asks for lat=256 chunks against a store written with lat=512, every read straddles two stored chunks and decompresses both. This gate asserts that on-disk encoding chunks divide the requested Dask chunks cleanly.

def assert_chunk_alignment(ds: xr.Dataset, requested: dict[str, int]) -> None:
    """
    Raise if requested Dask chunks are not integer multiples of the stored
    encoding chunks. Misalignment silently amplifies object-store reads.
    """
    for var in ds.data_vars:
        stored = ds[var].encoding.get("chunks")
        if stored is None:
            continue
        for dim, stored_len in zip(ds[var].dims, stored):
            if dim in requested and requested[dim] % stored_len != 0:
                raise ValueError(
                    f"'{var}' dim '{dim}': requested chunk {requested[dim]} is not a "
                    f"multiple of stored chunk {stored_len}; reads will amplify."
                )
    logger.info("Chunk alignment gate passed.")

Gate 3 — Consolidated Metadata Present #

A Zarr store without .zmetadata forces the reader to list and GET one metadata object per array — a “metadata storm” that can add seconds to job startup and dollars to the object-store bill at scale. Assert the consolidated file exists before publishing.

import fsspec


def assert_consolidated(store_uri: str, storage_options: dict | None = None) -> None:
    """Raise if a Zarr store lacks consolidated metadata."""
    fs, _, (root,) = fsspec.get_fs_token_paths(store_uri, storage_options=storage_options or {})
    if not fs.exists(f"{root.rstrip('/')}/.zmetadata"):
        raise RuntimeError(
            f"{store_uri} has no .zmetadata; re-run to_zarr(..., consolidated=True) "
            "or zarr.consolidate_metadata(store)."
        )
    logger.info("Consolidated-metadata gate passed: %s", store_uri)

Common Failure Modes and Diagnosis #

1. NetCDF-on-Object-Store Whole-File GETs #

Symptom: Opening a single NetCDF on S3 pulls hundreds of megabytes and stalls for seconds before any data is read; egress costs spike; a Dask cluster’s network saturates on job startup.

Root cause: A NetCDF4 file is one HDF5 container. The library must walk the internal B-tree to locate chunks, and a naive xr.open_dataset("s3://…nc") without a byte-range-capable driver downloads the entire file to a temp path first. There is no per-chunk object to fetch — the addressing granularity is the whole file.

Diagnosis: Watch object-store access logs or s3fs debug logging; a single open producing one large sequential GET (rather than many small range reads) confirms whole-file transfer. Fix by serving NetCDF through THREDDS/OPeNDAP, wrapping it with kerchunk to expose virtual chunk references, or converting to Zarr for the compute tier.

2. Un-consolidated Zarr Metadata Storms #

Symptom: xr.open_zarr takes 5-30 seconds before the first chunk is read; the object store shows hundreds of tiny GETs and LISTs at open time; latency scales with the number of variables.

Root cause: Without .zmetadata, the reader fetches .zgroup, .zattrs, and .zarray for every array and group individually. On high-latency object stores each of those is a separate round trip.

Diagnosis: Count GETs during open. If they scale with variable count, consolidated metadata is missing. Remediate:

import zarr

zarr.consolidate_metadata(store)          # rebuild .zmetadata in place
# then always read with:
# xr.open_zarr(store, consolidated=True)

The dedicated guide on fixing Zarr consolidated metadata on GCS covers the GCS-specific eventual-consistency wrinkles.

3. Chunk-Size Mismatch (Read Amplification) #

Symptom: A time-series extraction at one grid point reads gigabytes; Dask memory spikes far beyond the array footprint.

Root cause: The store is chunked for spatial access (time=1, lat=512, lon=512) but the workload extracts long time series at fixed (lat, lon). Every timestep lives in a different chunk, so a single-pixel time series touches every chunk in the store.

Diagnosis: Compare the workload’s access pattern to the stored chunk grid. If the iterated dimension is chunked at length 1, rechunk to a time-contiguous layout (time=8760, lat=64, lon=64) with rechunker before the analysis stage. This is the array analogue of the HDF5 chunk-alignment failure documented in the NetCDF vs GeoTIFF guide.

4. HDF5 SWMR / File Locking on Shared Storage #

Symptom: Concurrent NetCDF writers on an NFS or Lustre mount raise HDF5 error: unable to lock file or produce silently truncated files; readers see NetCDF: HDF error.

Root cause: HDF5 uses POSIX file locks that many networked filesystems implement incompletely. Single-writer/multiple-reader (SWMR) mode helps for one writer but does not enable the concurrent multi-writer pattern coastal fan-out jobs need. Object storage has no POSIX locks at all.

Diagnosis: If a parallel job intermittently corrupts a NetCDF output, it is almost always a locking or concurrent-write violation. The structural fix is to stop fighting HDF5’s single-writer model: write Zarr, where each worker owns disjoint chunk objects and no lock is required. Set HDF5_USE_FILE_LOCKING=FALSE only as a diagnostic, never as a production remedy — it masks the corruption rather than preventing it.

Pipeline Integration and Downstream Handoff #

Once a store passes all three gates, it is handed to two downstream stages. The Publishing COG & Zarr to Object Storage workflow takes the validated Zarr store URI, applies bucket lifecycle and cache-control policy, and registers it in the catalog; it also emits the manifest that records chunk grid, codec, and consolidated-metadata status so consumers do not have to introspect the store to configure their reads.

Large fan-out writes — a 40-year hindcast rechunked across many workers, or an ensemble surge forecast where each member writes its own region — are driven by Orchestrating Marine Pipelines with Dask & Prefect. That stage is where Zarr’s parallel region writes pay off: a Prefect flow pre-allocates the store’s time axis, then maps a task per region so every worker calls to_zarr(region=…) into disjoint chunks with no coordination beyond the pre-allocation. NetCDF cannot participate in that pattern, which is the single strongest reason cloud-native coastal pipelines standardise on Zarr internally and reserve NetCDF for the ERDDAP/THREDDS-facing archival boundary.

The manifest each store carries mirrors the schema used across this domain: it names the source, the CRS WKT, the chunk grid, the codec, the CF version, and the gates passed. Downstream stages read the manifest and abort if the store was published without passing the alignment gate — the same fail-fast contract that governs the CRS alignment handoff.

---

Related #

• Understanding NetCDF vs GeoTIFF for Marine Data — the raster fork: when a 2D grid belongs in a GeoTIFF/COG rather than an N-D array
• CRS Alignment for Coastal GIS Projects — the upstream stage that assigns the CRS and grid_mapping both writers assume
• Publishing COG & Zarr to Object Storage — bucket policy, cache-control, and catalog registration for validated stores
• Orchestrating Marine Pipelines with Dask & Prefect — fan-out region writes that exploit Zarr’s concurrent-write model

Up: Marine Spatial Data Fundamentals & Architecture