Lazy evaluation & Dask, explained¶
This is the concept notebook for pyramids' lazy path — the mental model, taught the way the
xarray + Dask tutorial teaches it. Once the model
clicks, the two reference cookbooks (lazy-dataset-complete.ipynb, lazy-collection-complete.ipynb) are just the
full API surface applied.
Why bother with Dask? It buys you two things:
- Parallelism — split the work across all your CPU cores (or a whole cluster) instead of one.
- Bigger-than-RAM data — process a raster (or a stack of them) that would never fit in memory, by streaming it one chunk at a time. Dask calls this out-of-core: you scale from memory-sized data to disk-sized data.
The one idea underneath both: a lazy array is a recipe, not the food. It records how to read and compute —
tile by tile — and cooks only when you ask. In pyramids the single switch that hands you a recipe instead of a
NumPy array is chunks=.
Caveat up front, straight from the xarray tutorial: using Dask does not always make your computation faster. Section 6 shows exactly when it pays off and when it just adds overhead.
Setup¶
We import numpy, dask.array, and a small helper that draws a dask array's chunk boundaries as an image (so the
tiling is visible even when your notebook is untrusted and hides dask's HTML repr).
%matplotlib inline
import time
from pathlib import Path
import numpy as np
import dask.array as da
import matplotlib.pyplot as plt
from pyramids.dataset import Dataset
from pyramids.base.protocols import as_numpy, is_lazy
DATA = (Path('..') / '..' / '..' / 'examples' / 'data').resolve()
def plot_chunk_grid(arr, title):
# Draw the dask chunk boundaries over the array's row/column extent.
ny, nx = arr.shape[-2], arr.shape[-1]
ys, xs = np.cumsum((0,) + arr.chunks[-2]), np.cumsum((0,) + arr.chunks[-1])
fig, ax = plt.subplots(figsize=(4.5, 4))
for y0, y1 in zip(ys[:-1], ys[1:]):
for x0, x1 in zip(xs[:-1], xs[1:]):
ax.add_patch(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, fill=False, ec='crimson', lw=1.4))
ax.set(xlim=(0, nx), ylim=(ny, 0), title=title, xlabel='columns', ylabel='rows')
ax.set_aspect('equal')
return fig
DATA.is_dir()
2026-07-11 14:38:32 | INFO | pyramids.base.config | Logging is configured.
True
1. Eager vs lazy — one switch¶
read_array() with no arguments is eager: it reads the whole band into a numpy.ndarray right now.
Add chunks= and the very same method returns a dask.array instead — the recipe. Nothing is read on that
call. is_lazy (pyramids' backend-agnostic check) tells the two apart.
ds = Dataset.read_file(DATA / 'geotiff' / 'south-america-mswep_1979010100.tif')
eager = ds.read_array() # numpy — read happens now
lazy = ds.read_array(chunks=(256, 256)) # dask — a recipe, nothing read yet
print("eager :", type(eager).__name__, eager.shape, "| is_lazy =", is_lazy(eager))
print("lazy :", type(lazy).__name__, lazy.shape, "| is_lazy =", is_lazy(lazy))
eager : ndarray (780, 850) | is_lazy = False lazy : Array (780, 850) | is_lazy = True
2. Chunks — the blocked algorithm¶
Dask turns one big array into a grid of small blocks (chunks) and runs every operation block-by-block. That is what lets an out-of-core array work: only one chunk needs to be in memory at a time. Displaying the lazy array shows dask's summary — array vs chunk size, the count, and an SVG of the layout:
lazy
|
||||||||||||||||
The same tiling as plain text and as an image (both render regardless of notebook trust). Our 780 × 850 band cut at 256 × 256 becomes a 4 × 4 grid of 16 chunks — the right column and bottom row are smaller remainders.
print("chunk grid :", lazy.numblocks, "->", lazy.npartitions, "chunks")
print("chunk size :", lazy.chunksize)
print("per-axis :", lazy.chunks)
chunk grid : (4, 4) -> 16 chunks chunk size : (256, 256) per-axis : ((256, 256, 256, 12), (256, 256, 256, 82))
plot_chunk_grid(lazy, f"dask chunk grid — {lazy.numblocks[-2]}x{lazy.numblocks[-1]} = {lazy.npartitions} tiles")
3. Lazy computation — building a graph¶
Here is the crux. Every operation you apply to a lazy array is deferred: instead of running, dask appends a
task to a graph that records the work to do. No pixel is read. We chain three operations — mask the -9999
fill, scale, then subtract the scene mean — and watch the task count grow while zero bytes are read:
nodata = float(np.ravel(ds.no_data_value)[0])
# Build the pipeline step by step — each stage extends the graph, none runs.
read = lazy
masked = da.where(lazy == nodata, np.nan, lazy)
scaled = masked * 2
anomaly = scaled - da.nanmean(scaled)
stages = {'1. read': read, '2. mask fill': masked, '3. scale x2': scaled, '4. anomaly': anomaly}
graph_sizes = {name: len(arr.__dask_graph__()) for name, arr in stages.items()}
graph_sizes # tasks queued after each step — still nothing computed
{'1. read': 16, '2. mask fill': 48, '3. scale x2': 64, '4. anomaly': 101}
Each step adds tasks rather than doing work. Plotting the growth makes the point — the graph inflates from a handful of reads to dozens of queued tasks, and not a single one has run yet:
fig, ax = plt.subplots(figsize=(5.5, 3))
ax.bar(list(graph_sizes.keys()), list(graph_sizes.values()), color='steelblue')
ax.set(ylabel='tasks in graph', title='the graph grows as you compose — nothing has run')
for i, v in enumerate(graph_sizes.values()):
ax.text(i, v + 1, str(v), ha='center')
fig.autofmt_xdate(rotation=20)
plt.show()
The graph is organised into named layers, one per operation — dask's record of the pipeline it will execute:
[name.split('-')[0] for name in anomaly.dask.layers] # the op behind each layer
['sub', 'where', '_read_chunk', 'eq', 'mul', 'mean_chunk', 'mean_combine', 'mean_agg']
4. Getting concrete values — compute, as_numpy, persist¶
The graph runs only when you ask. There are three ways to ask, and the difference matters (they mirror xarray's
.compute() / .load() / .persist()):
| Call | Runs the graph? | Returns | Keeps chunking? |
|---|---|---|---|
arr.compute() |
yes | a numpy array | no — fully materialised; arr itself is untouched and reusable |
as_numpy(arr) |
yes | a numpy array | no — pyramids' backend-agnostic "give me numpy" (no-op on numpy input) |
arr.persist() |
yes | a dask array, results cached in RAM | yes — stays lazy-shaped, but re-using it won't recompute |
computed = anomaly.compute() # -> numpy, the anomaly array is unchanged and can be reused
via_helper = as_numpy(anomaly) # -> numpy, same result via the pyramids helper
persisted = anomaly.persist() # -> still a dask array, but its chunks are now in memory
print("compute() ->", type(computed).__name__, computed.shape)
print("as_numpy ->", type(via_helper).__name__, via_helper.shape)
print("persist() ->", type(persisted).__name__, "(still dask, chunks cached)")
print("anomaly is still lazy afterwards:", is_lazy(anomaly))
compute() -> ndarray (780, 850) as_numpy -> ndarray (780, 850) persist() -> Array (still dask, chunks cached) anomaly is still lazy afterwards: True
5. Watch it run — a progress bar¶
During section 3 nothing happened; compute() is where the 16 tile-reads and the reduction actually stream through
the scheduler. dask.diagnostics.ProgressBar makes that visible — the bar fills as chunks complete:
from dask.diagnostics import ProgressBar
with ProgressBar():
result = anomaly.compute()
float(np.nanmin(result)), float(np.nanmax(result))
[ ] | 0% Completed | 187.33 us
[########################################] | 100% Completed | 102.55 ms
(-8.554137229919434, 274.6658630371094)
6. Does Dask always help?¶
No — and pretending otherwise sets you up to be disappointed. For a small raster that already fits in RAM, the scheduler's bookkeeping usually makes the lazy path slower than a plain NumPy read. We time both on this 660 k-pixel band:
def eager_anomaly():
a = ds.read_array().astype('float32')
a = np.where(a == nodata, np.nan, a) * 2
return a - np.nanmean(a)
def lazy_anomaly():
a = ds.read_array(chunks=(256, 256))
a = da.where(a == nodata, np.nan, a) * 2
return (a - da.nanmean(a)).compute()
def _best(fn, n=3):
best = float('inf')
for _ in range(n):
t0 = time.perf_counter()
fn()
best = min(best, time.perf_counter() - t0)
return best
print(f"eager (numpy) : {_best(eager_anomaly) * 1e3:6.1f} ms")
print(f"lazy (dask) : {_best(lazy_anomaly) * 1e3:6.1f} ms")
print("On data this small, dask's overhead usually loses — that's expected.")
eager (numpy) : 6.6 ms
lazy (dask) : 28.7 ms On data this small, dask's overhead usually loses — that's expected.
So when does the lazy path pay off? Reach for chunks= when:
- the raster (or the stack) is bigger than RAM — chunking is the only way to touch it at all;
- the work is CPU-heavy and parallel (focal filters, DEM derivatives, per-tile band math) and you have cores to spread it over;
- the data lives on the cloud (COGs on S3) — dask fetches only the byte-ranges of the chunks it needs;
- you are reducing a multi-file time cube and want each timestep read on a different worker.
Skip it for small rasters that comfortably fit in memory: a plain read_array() is simpler and faster.
Further reading¶
- xarray + Dask tutorial — the model this notebook mirrors.
- Dask best practices — chunk sizing and the pitfalls.
lazy-dataset-complete.ipynb— every lazyDatasetsurface (focal ops, DEM derivatives, Zarr I/O, zonal stats).lazy-collection-complete.ipynb— the 4-DDatasetCollectioncube (reductions, groupby, parallel Zarr write).dask-lazy-datasets.ipynb— lazy reads +map_blocks+ parallelDelayedwrites end to end.