Building a YAML-Driven Multi-Index Pipeline

Instead of hard-coding NDVI, EVI, SAVI, and NDWI as separate functions, declare each one in a YAML registry that maps a formula string to band tokens, then evaluate every formula through a single safe engine over xarray DataArrays:

# indices.yaml
ndvi:
  formula: "(nir - red) / (nir + red)"
  bands: [nir, red]
ndwi:
  formula: "(green - nir) / (green + nir)"
  bands: [green, nir]
savi:
  formula: "1.5 * (nir - red) / (nir + red + 0.5)"
  bands: [nir, red]

This page belongs to Spectral Index Calculation Pipelines and focuses on the registry-plus-engine pattern. The underlying array arithmetic is exactly the vectorised Band Math Operations with xarray; what is new here is decoupling the formulas from the code so a scientist can add an index by editing YAML rather than shipping a release.

Why This Arises in Remote Sensing Workflows

A team that starts with one NDVI function soon needs EVI, then SAVI for arid sites, then NDWI for water, then a half-dozen tasselled-cap and red-edge indices. Copy-pasting a function per index produces near-duplicate code, drifting band conventions, and a codebase where adding an index means a code review and a deploy. It also couples the index definitions to one sensor’s band naming.

A registry inverts that. The formulas live as data — a YAML file listing each index’s expression and the band tokens it consumes — and a single engine loads the scene’s bands, binds the tokens, and evaluates every formula. Adding an index becomes a one-line YAML edit that a domain scientist can make without opening the codebase, and the set of supported indices is now something you can diff, review, and version independently of the processing engine.

The only real engineering problem is evaluating a formula string safely. Python’s built-in eval would run the string as code, so a formula like __import__('os').system('...') slipped into a YAML file — by an attacker or simply by a typo that happens to be valid Python — would execute with the pipeline’s privileges. The engine sidesteps this entirely by parsing the expression with the ast module into a syntax tree and then walking only the nodes it explicitly permits: numeric constants, names that resolve to known bands, and the four or five arithmetic operators. Anything else — a function call, an attribute access, a subscript, a comparison — raises before any evaluation happens. The result is an evaluator with the ergonomics of eval but a tightly bounded blast radius.

This topic sits in Spectral Index Calculation Pipelines within Satellite Processing Workflows & Index Pipelines. It pairs naturally with masking — feed it the clean reflectance from Masking Clouds with the Sentinel-2 SCL Band so every index is computed only over valid pixels.

How the Registry and Engine Fit Together

The diagram traces one scene through the registry: bands load into a Dataset, each formula is parsed to a validated expression tree, and the engine evaluates them into a stack of index layers.

Because the engine binds band tokens per formula, indices with different band requirements — NDVI on red/nir, NDWI on green/nir — flow through the same evaluator with no special-casing. The Dataset of loaded bands acts as the shared namespace: each formula names the subset it needs, the engine checks those names exist, and evaluation proceeds only when the dependency is satisfied. That check is what turns a vague runtime failure into a clear, up-front error naming the missing band and the index that wanted it.

The output is a Dataset with one variable per index, all sharing the scene’s spatial coordinates, so it serialises cleanly to NetCDF or Zarr and slots straight into any downstream step that expects labelled geospatial arrays. Nothing about the engine is Sentinel-2 specific; point the band tokens at Landsat, PlanetScope, or a drone orthomosaic and the same registry produces the same indices.

Environment & Setup

pip install "xarray>=2023.6" "rioxarray>=0.15" "PyYAML>=6.0" "numpy>=1.24"

Complete Working Example

The engine parses each formula into an AST, validates that it contains only numbers, known band names, and arithmetic operators, then evaluates the tree by binding band tokens to DataArrays.

import ast
import operator
from pathlib import Path
import yaml
import xarray as xr
import rioxarray  # registers the .rio accessor

# Only these node/operator types are allowed in a formula
_BIN_OPS = {
    ast.Add: operator.add, ast.Sub: operator.sub,
    ast.Mult: operator.mul, ast.Div: operator.truediv,
    ast.Pow: operator.pow,
}
_UNARY_OPS = {ast.UAdd: operator.pos, ast.USub: operator.neg}

def _eval_node(node: ast.AST, bands: dict[str, xr.DataArray]) -> xr.DataArray:
    """Recursively evaluate a validated arithmetic AST over band DataArrays."""
    if isinstance(node, ast.Expression):
        return _eval_node(node.body, bands)
    if isinstance(node, ast.Constant):          # a numeric literal, e.g. 0.5
        if not isinstance(node.value, (int, float)):
            raise ValueError("Only numeric constants are allowed.")
        return node.value
    if isinstance(node, ast.Name):              # a band token, e.g. nir
        if node.id not in bands:
            raise KeyError(f"Unknown band '{node.id}' in formula.")
        return bands[node.id]
    if isinstance(node, ast.BinOp) and type(node.op) in _BIN_OPS:
        return _BIN_OPS[type(node.op)](
            _eval_node(node.left, bands), _eval_node(node.right, bands)
        )
    if isinstance(node, ast.UnaryOp) and type(node.op) in _UNARY_OPS:
        return _UNARY_OPS[type(node.op)](_eval_node(node.operand, bands))
    raise ValueError(f"Disallowed expression element: {ast.dump(node)}")

def evaluate_formula(formula: str, bands: dict[str, xr.DataArray]) -> xr.DataArray:
    """Safely evaluate one index formula string over the given bands."""
    tree = ast.parse(formula, mode="eval")   # never exec; parse only
    return _eval_node(tree, bands)

def run_index_pipeline(
    registry_path: str,
    band_paths: dict[str, str],
) -> xr.Dataset:
    """
    Compute every index in the YAML registry from the scene's bands.
    band_paths maps a band token (e.g. 'nir') to a GeoTIFF path.
    """
    registry = yaml.safe_load(Path(registry_path).read_text())

    # Load each band once as a squeezed 2-D DataArray
    bands = {
        token: rioxarray.open_rasterio(path).squeeze("band", drop=True).astype("float32")
        for token, path in band_paths.items()
    }

    out = xr.Dataset()
    for name, spec in registry.items():
        needed = spec["bands"]
        missing = [b for b in needed if b not in bands]
        if missing:
            raise KeyError(f"Index '{name}' needs missing bands: {missing}")
        out[name] = evaluate_formula(spec["formula"], bands)
    return out


if __name__ == "__main__":
    indices = run_index_pipeline(
        registry_path="indices.yaml",
        band_paths={
            "red": "B04_10m.tif",
            "green": "B03_10m.tif",
            "nir": "B08_10m.tif",
        },
    )
    print(list(indices.data_vars))          # ['ndvi', 'ndwi', 'savi', ...]
    indices.to_netcdf("scene_indices.nc")

The critical safety property is that ast.parse(..., mode="eval") builds a tree without executing anything, and _eval_node only descends into node types on its allow-list. A formula containing a function call, an attribute access, or an import raises ValueError instead of running.

Variant Patterns

1. Per-sensor band aliasing

Keep formulas sensor-agnostic by resolving band tokens through an alias map, so the same indices.yaml runs on Sentinel-2 and Landsat.

SENTINEL2 = {"red": "B04", "green": "B03", "nir": "B08", "swir1": "B11"}
LANDSAT9  = {"red": "SR_B4", "green": "SR_B3", "nir": "SR_B5", "swir1": "SR_B6"}

def resolve(band_files: dict[str, str], alias: dict[str, str]) -> dict[str, str]:
    return {token: band_files[code] for token, code in alias.items() if code in band_files}

2. Declaring constants and coefficients in the registry

EVI needs gain and soil coefficients. Put them in the formula directly, or extend each entry with a constants block the engine binds alongside the bands.

evi:
  formula: "2.5 * (nir - red) / (nir + 6*red - 7.5*blue + 1)"
  bands: [nir, red, blue]

3. Streaming over a Dask-backed Dataset

Open the bands with chunks= so the identical engine evaluates lazily and scales to full tiles. The arithmetic tree is unchanged; xarray builds a task graph instead of computing eagerly.

import rioxarray
bands = {
    t: rioxarray.open_rasterio(p, chunks={"x": 2048, "y": 2048}).squeeze("band", drop=True)
    for t, p in band_paths.items()
}
# evaluate_formula(...) now returns a lazy DataArray; call .compute() to run

Common Errors

ValueError: Disallowed expression element

The formula contains something outside plain arithmetic — a function call like sqrt(...), an attribute, or a comparison. Restrict formulas to +, -, *, /, **, parentheses, band names, and numeric constants, or extend the allow-list deliberately.

KeyError: Unknown band 'swir2' in formula

A formula references a band token that was not loaded into the Dataset. Confirm every token in the registry’s bands list is present in band_paths (or resolvable through your alias map) before running.

Indices come out all NaN or all zero

The bands were read as integers and an integer division truncated, or a scale factor was not applied. Cast bands to float32 on load and apply any reflectance scale before evaluating the formulas.

Related

• Spectral Index Calculation Pipelines — the parent section on structuring multi-index computation.
• Band Math Operations with xarray — the vectorised DataArray arithmetic the engine relies on.
• Masking Clouds with the Sentinel-2 SCL Band — mask reflectance before feeding it to the index engine.
• Computing EVI and NDWI from Sentinel-2 Bands — the explicit per-index formulas the registry generalises.