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

What opening reads, what stays lazy, what forces a read

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.

Lazy against loaded on a Snowdon Sentinel-2 cube

The data used on this page

cube.zarr is built from every Sentinel-2 L2A scene in Earth Search over a 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N) from 1 June to 31 December 2025 — 127 scenes, the same stack as How to build a cloud-free composite in Python. It holds reflectance (time, band, y, x) with the red and near-infrared bands, stored as integers with a 1/10,000 scale factor, and the scene classification scl, one scene per chunk.

import os

import numpy as np
import pystac_client
import rasterio
import xarray as xr
from pyproj import Transformer
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
if not os.path.exists("cube.zarr"):
    to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
    bbox = [*to_wgs84.transform(*BOUNDS[:2]), *to_wgs84.transform(*BOUNDS[2:])]
    items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
        collections=["sentinel-2-l2a"], bbox=bbox, datetime="2025-06-01/2025-12-31").items(),
        key=lambda i: i.id)

    def window(href):
        with rasterio.open(href) as src:
            return src.read(1, window=from_bounds(*BOUNDS, transform=src.transform))

    os.makedirs("scenes", exist_ok=True)
    for n, item in enumerate(items):
        path = f"scenes/{item.id}.npz"
        if not os.path.exists(path):
            scl = window(item.assets["scl"].href).repeat(2, axis=0).repeat(2, axis=1)
            np.savez_compressed(path, red=window(item.assets["red"].href), nir=window(item.assets["nir"].href),
                                scl=scl, date=str(item.datetime.date()))
        z = np.load(path)
        ds = xr.Dataset(
            {"reflectance": (("time", "band", "y", "x"), np.stack([z["red"], z["nir"]])[None] / np.float32(10_000)),
             "scl": (("time", "y", "x"), z["scl"][None].astype("uint8"))},
            coords={"time": [np.datetime64(str(z["date"]))], "band": ["red", "nir"],
                    "y": 5_885_435 - 10 * np.arange(1000), "x": 422_885 + 10 * np.arange(1000)})
        ds.reflectance.encoding = {"scale_factor": np.float32(1e-4), "dtype": "uint16", "_FillValue": 0}
        ds.to_zarr("cube.zarr", mode="w" if n == 0 else "a", zarr_format=2,
                  **({} if n == 0 else {"append_dim": "time"}))
print(xr.open_zarr("cube.zarr"))

<xarray.Dataset> Size: 2GB
Dimensions:      (time: 127, band: 2, y: 1000, x: 1000)
Coordinates:
  * time         (time) datetime64[ns] 1kB 2025-06-04 2025-06-07 ... 2025-12-29
  * band         (band) <U3 24B 'red' 'nir'
  * y            (y) int64 8kB 5885435 5885425 5885415 ... 5875455 5875445
  * x            (x) int64 8kB 422885 422895 422905 ... 432855 432865 432875
Data variables:
    reflectance  (time, band, y, x) float64 2GB dask.array<chunksize=(1, 1, 500, 500), meta=np.ndarray>
    scl          (time, y, x) uint8 127MB dask.array<chunksize=(1, 500, 500), meta=np.ndarray>

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


cube = xr.open_zarr("cube.zarr")
inspect_laziness(cube["reflectance"], "reflectance")

  reflectance: lazy (dask)
    shape (127, 2, 1000, 1000), dtype float64
    2.03 GB if fully loaded
    chunks (1, 1, 500, 500), 1016 partitions
    1,017 tasks in the graph

Printing the fully loaded size before computing is the habit that prevents most out-of-memory failures. Two gigabytes is fine on a workstation; the same code on a continental cube is not. The integers on disk decode to float64 through the scale factor, which is why the loaded size is four times the stored one.

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


bands = cube.assign(red=cube.reflectance.sel(band="red"), nir=cube.reflectance.sel(band="nir"))
composite = masked_composite(bands)

  graph has 9,165 tasks before compute
  median NDVI 0.612, observations 4-31

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)


nir = cube.reflectance.sel(band="nir").chunk({"time": -1, "y": 250, "x": 250})
lazy = apply_chunkwise(nir, p90)
print(f"  computed: NIR 90th percentile over time, median {float(lazy.compute().median()):.3f}")

  lazy result: (1000, 1000), chunks (250, 250)
  computed: NIR 90th percentile over time, median 1.003

The median of 1.003 is above 1 because nothing here was masked: the brightest tenth of 127 scenes over Snowdon is mostly cloud. 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.

Why laziness is a necessity, and its traps

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, 1, 500, 500) — 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.

Forced reads, .values in loops and persist

• 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.

Internal links

• Chunked arrays and Zarr explained — where the chunks come from
• How to turn a STAC search into an xarray data cube — building a lazy cube from imagery
• How to choose chunk and tile sizes that actually help — sizing them
• Cloud-native geospatial explained — why partial reads are possible at all
• How to scale a GeoPandas job with dask-geopandas — the same laziness for tables
• dask-geopandas uses more memory or time than plain GeoPandas — when laziness costs more than it saves
• How to process a very large GeoPackage in chunks with Python — manual chunking
• How to profile a slow Python GIS script and find the real bottleneck — finding the accidental compute
• How to open a NetCDF file in Python with xarray — the lazy open in practice
• How to open hundreds of NetCDF files as one dataset — laziness across many files

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.