Computing Band Ratios with Dask-Backed DataArrays

Open each band with identical chunks, write the ratio as ordinary arithmetic, and let the write drive the computation:

import rioxarray

chunks = {"x": 2048, "y": 2048}
red = rioxarray.open_rasterio("B04.tif", chunks=chunks, mask_and_scale=True).squeeze("band")
nir = rioxarray.open_rasterio("B08.tif", chunks=chunks, mask_and_scale=True).squeeze("band")

ndvi = (nir - red) / (nir + red)              # a graph, not an array
ndvi.rio.to_raster("ndvi.tif", tiled=True, lock=False, compress="zstd")

Nothing is read until the write starts, and at no point does the whole scene sit in memory. This page belongs to band math operations with xarray in Core Raster Fundamentals & STAC Mapping.

What Dask Is Doing Underneath

Four chunks, four independent chains Each chunk of the output depends only on the matching chunk of the red band and the matching chunk of the near-infrared band. The graph is therefore four independent chains of read, read, subtract, add, divide and write, which Dask runs in parallel with peak memory of a few chunks rather than the whole scene. An embarrassingly parallel graph read red chunk 0 read nir chunk 0 (nir - red)/(nir + red) write chunk 0 chunk 1: the same chain, in parallel chunk 2: the same chain, in parallel chunk 3: the same chain, in parallel peak memory ≈ workers × 3 chunks not the scene

Because no chunk depends on any other, the graph is as parallel as it can be and the memory high-water mark is set by the number of chunks in flight, not by the scene size. That is the property that lets the same three lines run on a laptop for one tile and on a cluster for a continent.

Environment & Setup

Package Version pin Used for
rioxarray >=0.15 Chunked, masked-and-scaled reads and lock-free writes
xarray >=2024.1 Lazy labelled arithmetic
dask[array] >=2024.1 The chunked execution engine
distributed >=2024.1 Optional: a dashboard and a local cluster 
pip install "rioxarray>=0.15" "xarray>=2024.1" "dask[array]>=2024.1" "distributed>=2024.1"

Complete Working Example

import rioxarray
import xarray as xr
from dask.distributed import Client, LocalCluster


def open_band(path: str, chunk: int = 2048) -> xr.DataArray:
    da = rioxarray.open_rasterio(path, chunks={"x": chunk, "y": chunk},
                                 mask_and_scale=True, lock=False)
    return da.squeeze("band", drop=True)


def ratio_to_cog(a_path: str, b_path: str, out_path: str, *,
                 chunk: int = 2048) -> None:
    """(a - b) / (a + b) for two single-band rasters, streamed to a tiled GeoTIFF."""
    a, b = open_band(a_path, chunk), open_band(b_path, chunk)

    if a.rio.transform() != b.rio.transform() or a.shape != b.shape:
        b = b.rio.reproject_match(a)                   # align grids first
    b = b.chunk(a.chunks)                              # identical chunk boundaries

    den = a + b
    ratio = ((a - b) / den).where(den != 0)            # zero denominator -> NaN
    ratio = ratio.astype("float32").rio.write_nodata(float("nan"))

    ratio.rio.to_raster(
        out_path, tiled=True, blockxsize=512, blockysize=512,
        compress="zstd", predictor=3, lock=False,       # lock=False: parallel writes
        BIGTIFF="IF_SAFER",
    )


if __name__ == "__main__":
    with LocalCluster(n_workers=4, threads_per_worker=2, memory_limit="4GB") as cluster, \
         Client(cluster) as client:
        print("dashboard:", client.dashboard_link)
        ratio_to_cog("S2A_36NYF_B08.tif", "S2A_36NYF_B04.tif", "ndvi.tif")

lock=False on both read and write matters for throughput: by default rioxarray serialises GDAL access through a lock, which is safe and slow. With independent chunks and a tiled output, parallel access is safe, and removing the lock lets every worker read and write at once. The mechanics of the cluster itself are covered in scaling raster processing with Dask.

Chunk Alignment Across Bands

Misaligned chunks force a shuffle When both bands use 2048 pixel chunks, output chunk i needs only input chunk i from each band. When one band uses 2048 and the other 1830 — for example a 20 metre band resampled without matching chunks — every output chunk needs pieces of two input chunks from the second band, and Dask inserts a rechunk step with extra memory and network traffic. aligned chunks misaligned chunks boundaries line up; no shuffle every output chunk straddles two inputs Rechunk one band to the other's chunks explicitly, once, before the arithmetic.

Misalignment usually arrives through a resampling step: a 20 m band upsampled to 10 m keeps its own chunking, which no longer matches the native 10 m band. Dask handles it correctly, by inserting a rechunk that moves data between workers, but the cost shows up as memory spikes and slow progress in the dashboard. Calling b.chunk(a.chunks) once, explicitly, makes the cost visible and pays it only once. The deeper treatment of chunk choice is in tuning Dask chunk sizes for raster cubes.

Computing Several Ratios at Once

When a pipeline needs NDVI, NDWI and NDMI from the same bands, computing them separately reads every band three times. Building all three into one graph and writing them together lets Dask share the reads:

import dask
import xarray as xr

b03, b04, b08, b11 = (open_band(p) for p in ("B03.tif", "B04.tif", "B08.tif", "B11.tif"))
b11 = b11.rio.reproject_match(b08).chunk(b08.chunks)

def nd(x, y):
    den = x + y
    return ((x - y) / den).where(den != 0).astype("float32")

stack = xr.Dataset({"ndvi": nd(b08, b04), "ndwi": nd(b03, b08), "ndmi": nd(b08, b11)})
dask.compute(*[stack[v].rio.to_raster(f"{v}.tif", tiled=True, lock=False,
                                      compress="zstd", compute=False)
               for v in stack.data_vars])

Passing compute=False to each write returns a delayed object, and computing them together in one call means each band chunk is read once and consumed by all three indices. On remote data, where reads dominate, that alone can make a multi-index job two to three times faster.

Verification

Confirm lazy, aligned, correct Before writing, confirm the result is Dask-backed so nothing was loaded and that its chunks match the inputs. After writing, compare one window against an eager NumPy computation of the same ratio; they should match to floating-point precision. Three quick checks still lazy data is a Dask array chunks aligned same as the inputs matches eager one window, allclose The eager comparison is the one that catches a mask or alignment mistake. 
import dask.array as da
import numpy as np
import rasterio

assert isinstance(ratio.data, da.Array), "result was computed eagerly"
assert ratio.chunks == a.chunks, "chunks drifted from the inputs"

with rasterio.open("B08.tif") as n, rasterio.open("B04.tif") as r, \
     rasterio.open("ndvi.tif") as out:
    win = ((1000, 1256), (1000, 1256))
    nv = n.read(1, window=win).astype("float32") * 1e-4
    rv = r.read(1, window=win).astype("float32") * 1e-4
    eager = (nv - rv) / (nv + rv)
    lazy = out.read(1, window=win)
assert np.allclose(eager, lazy, equal_nan=True, atol=1e-6)

The eager comparison on one window is cheap and catches the bugs that laziness hides: a mask applied in one path and not the other, or a grid mismatch that the lazy path resolved differently from the eager one.

Common Errors

Memory climbs until workers are killed

A .compute() or .values materialised the full result. Write directly with to_raster, which streams chunk by chunk.

The write is no faster than a single thread

The default GDAL lock serialised every access. Pass lock=False to reads and writes on tiled files.

Every task shows a rechunk in the dashboard

Band chunks differ. Rechunk one band to the other’s chunks once, before the arithmetic.

Output has seams at chunk boundaries

A neighbourhood operation crept into the graph without overlap. Pure ratios cannot seam; check for a focal step and use map_overlap for it.

Frequently Asked Questions

Q: Why must every band have the same chunks? Because element-wise arithmetic pairs chunk i of one band with chunk i of the other. If the chunk boundaries differ, Dask must rechunk one band to match before it can compute, which adds a shuffle step and extra memory to every operation.

Q: When should I call persist? Only when an intermediate is reused several times and fits in cluster memory. For a single ratio written straight to disk, persisting forces the whole result into memory for no benefit; let the write drive the computation instead.

Q: What chunk size suits a Sentinel-2 band? A multiple of the file’s internal tile size, typically 1024 or 2048 pixels square for 10 metre bands. That gives chunks of a few to a few tens of megabytes and means each chunk reads whole internal tiles rather than fragments.

Q: Is Dask worth it for a single scene on a laptop? Often yes, because it bounds memory rather than because it is faster. A scene that does not fit in RAM as float32 simply cannot be processed eagerly; chunked, it runs in a fixed memory budget at close to disk speed.