🧭 Concept Cloud-native formats Performance & scale Cloud-Native GIS · step 8 of 14 Updated

Lazy Loading Explained: Why xarray Reads Nothing Until You Ask

Problem statement

xr.open_zarr on a terabyte returns instantly. xr.open_dataset on a directory of a thousand GeoTIFFs prints a full description of a dataset it has not read.

Then one innocuous line — a .plot(), a .mean(), a print() of the values — takes twenty minutes and exhausts memory.

Lazy loading is why both happen. Understanding where the boundary lies between "building a description" and "actually reading" is the difference between a pipeline that scales and one that dies on the first real dataset.

Quick answer

Nothing is read until something forces it:

import xarray as xr

cube = xr.open_zarr("cube.zarr")          # metadata only, instant
subset = cube.sel(time="2025-07")         # still nothing read
mean = subset["reflectance"].mean("time") # still nothing read

values = mean.compute()                    # NOW it reads

Four operations that force a read, and are easy to trigger by accident:

array.values      # materialises the whole array
array.compute()   # explicit
float(array)      # scalar conversion
array.plot()      # needs numbers
np.asarray(array) # any numpy function
Operations building a task graph with no data read, then a single compute call executing the whole graph.
Everything before `.compute()` is bookkeeping. Everything at it is I/O.

Step-by-step solution

1. Know what opening actually does

open_zarr and open_dataset read the metadata: dimensions, coordinates, dtypes, chunk shapes and attributes. That is kilobytes regardless of the array size.

The printed representation shows shapes and dtypes because those are metadata. It shows no values, which is the visible sign that nothing was read.

2. Know which operations stay lazy

Indexing, selection, arithmetic, reductions and reshaping all build a graph rather than executing:

cube.sel(time=slice("2025-06", "2025-09"))
cube["red"] / cube["nir"]
cube.mean("time")
cube.where(cube["scl"] != 9)

Each returns a new lazy object describing what to do. The cost is a graph node, not a byte of I/O.

3. Know which operations force a read

Anything needing actual numbers:

  • .values, .to_numpy(), np.asarray()
  • .compute(), .load(), .persist()
  • float(), int(), bool() on a scalar
  • .plot(), .to_netcdf(), .to_zarr()
  • printing a value, though not printing the object

The dangerous one is .values inside a loop. Each call re-executes the whole graph from the beginning, so a loop of a hundred iterations reads the data a hundred times.

4. Compute once, at the end 

result = (cube.sel(time=window)
              .where(cube["scl"].isin(GOOD))
              .median("time"))
result.compute()                # one pass over the data

Building the whole expression before computing lets the scheduler read each chunk once and apply everything to it. Computing intermediate steps forces a separate pass for each.

5. Use persist for something reused 

masked = cube.where(mask).persist()      # compute and keep in memory
a = masked.mean("time").compute()
b = masked.max("time").compute()

persist executes the graph and keeps the result in memory as chunks. Without it, a and b each re-read and re-mask everything.

A loop calling .values on a lazy array re-executing the whole graph each iteration, against computing once outside the loop.
`.values` inside a loop is the most expensive line in most xarray code.

Code examples

Example 1 — telling lazy from loaded 

import xarray as xr
import numpy as np


def inspect_laziness(obj, name="array"):
    """Is this backed by dask, and how big would loading it be?"""
    data = obj.data if hasattr(obj, "data") else obj
    is_dask = hasattr(data, "dask")

    nbytes = obj.nbytes if hasattr(obj, "nbytes") else data.nbytes
    print(f"  {name}: {'lazy (dask)' if is_dask else 'loaded (numpy)'}")
    print(f"    shape {tuple(obj.shape)}, dtype {obj.dtype}")
    print(f"    {nbytes / 1e9:.2f} GB if fully loaded")

    if is_dask:
        print(f"    chunks {data.chunksize}, {data.npartitions} partitions")
        print(f"    {len(data.dask):,} tasks in the graph")
        if nbytes > 8e9:
            print("    ! do not call .values on this")
    return is_dask
  reflectance: lazy (dask)
    shape (120, 4, 509, 543), dtype float32
    0.53 GB if fully loaded
    chunks (1, 4, 509, 543), 120 partitions
    481 tasks in the graph

Printing the fully loaded size before computing is the habit that prevents most out-of-memory failures. Half a gigabyte is fine; the same code on a continental cube is not.

Example 2 — the pattern that reads once 

import xarray as xr
import numpy as np

UNUSABLE = [0, 1, 3, 8, 9, 10]


def masked_composite(cube, red="red", nir="nir", scl="scl"):
    """Build the whole expression, then compute it once."""
    usable = ~cube[scl].isin(UNUSABLE)

    red_masked = cube[red].where(usable)
    nir_masked = cube[nir].where(usable)
    ndvi = (nir_masked - red_masked) / (nir_masked + red_masked)

    result = xr.Dataset({
        "ndvi_median": ndvi.median("time"),
        "observations": usable.sum("time"),
    })

    print(f"  graph has {len(result['ndvi_median'].data.dask):,} tasks "
          "before compute")
    computed = result.compute()
    print(f"  median NDVI {float(computed['ndvi_median'].median()):.3f}, "
          f"observations {int(computed['observations'].min())}"
          f"-{int(computed['observations'].max())}")
    return computed

Both outputs are in one Dataset before .compute(), so the scheduler reads each chunk once and produces both. Computing them separately would double the I/O.

Example 3 — chunk-wise processing that never materialises 

import xarray as xr
import numpy as np


def apply_chunkwise(cube, func, output_dtype="float32"):
    """Run a NumPy function on each chunk without loading the whole array."""
    result = xr.apply_ufunc(
        func, cube,
        input_core_dims=[["time"]],
        output_dtypes=[output_dtype],
        dask="parallelized",
        vectorize=False,
    )
    print(f"  lazy result: {tuple(result.shape)}, "
          f"chunks {result.data.chunksize}")
    return result


def p90(values, axis=-1):
    return np.nanpercentile(values, 90, axis=axis)

dask="parallelized" is what keeps this lazy. Without it, apply_ufunc computes the input first, which defeats the purpose — and the failure mode is memory exhaustion rather than an error.

input_core_dims tells dask which dimension the function consumes, so it can rechunk to make that dimension contiguous before applying the function.

Explanation

Why laziness is necessary, not clever

A continental Sentinel-2 cube is petabytes. Any API that reads on open is unusable for it.

Laziness inverts the model: describe the computation first, then execute only what the result requires. Selecting one month from ten years reads one month's chunks, because the selection happened before any I/O.

That is the same argument as for lazy evaluation in databases, and it has the same consequence: the cost of an operation is not where it appears in the code.

Why .values in a loop is so expensive

A lazy object holds a graph, not a cache. Calling .values executes the graph and returns an array; calling it again executes the graph again.

Inside a loop, that means re-reading and re-computing everything on every iteration. A ten-line loop can therefore take a hundred times longer than the same work done once outside it, with no visible difference in the code.

The fix is always the same: compute once outside the loop, or .persist() if the result is reused several times in different ways.

Why chunk alignment matters

The chunk shape decides how much data each task reads. If the stored chunks are (1, 4, 509, 543) — one time step each — then a reduction over time must touch every chunk, which is correct and unavoidable.

But if the dask chunks are set differently from the stored chunks, each task reads parts of several stored chunks, and each stored chunk is read by several tasks. The same bytes are fetched more than once.

Passing chunks= matching the storage layout, or a whole multiple of it, avoids that.

Why the graph itself can be a problem

Dask graphs are Python objects, and a graph with millions of tasks costs real memory and scheduling time.

The usual cause is very small chunks: a cube with a million chunks and ten operations gives ten million tasks. The scheduler then spends more time managing the graph than reading data.

If a computation is slow with low CPU and low network use, check the task count. It is usually a chunking problem rather than a compute one.

A lazy pipeline building mask, index and two reductions before a single compute call reads the data once.
Two `.compute()` calls read everything twice. One `dask.compute(a, b)` reads it once.

Edge cases or notes

  • Opening reads metadata only, regardless of array size.
  • .values, .compute(), .plot() and float() all force a read.
  • .values inside a loop re-executes the whole graph every iteration.
  • Build the whole expression, then compute once.
  • .persist() for a result reused several ways.
  • dask="parallelized" in apply_ufunc, or it computes the input first.
  • Align dask chunks with storage chunks.
  • Check the task count when a computation is slow with idle CPU.

FAQ

Why does opening a huge dataset take no time?

Because only the metadata is read — dimensions, dtypes, chunk shapes. That is kilobytes regardless of the array size.

What forces xarray to actually read the data?

.values, .compute(), .load(), .plot(), float() on a scalar, writing to a file, and any NumPy function applied to the object.

Why is my loop over an xarray object so slow?

You are probably calling .values inside it. Each call re-executes the whole graph, so the data is read once per iteration.

What is the difference between compute and persist?

compute returns a concrete result. persist executes the graph and keeps the result in memory as chunks, so several later operations reuse it instead of re-reading.

Should my dask chunks match the storage chunks?

Yes, or be whole multiples of them. Misaligned chunks cause each stored chunk to be fetched more than once.

Why is my computation slow with idle CPU?

Often too many tasks. Very small chunks produce enormous graphs, and the scheduler spends its time on bookkeeping rather than data.

Do I need dask for xarray?

No. Without chunks=, arrays are NumPy-backed and load on access. Use dask when the data is larger than memory.