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.

YAML registry plus band Dataset feeding a safe evaluator A YAML registry of formulas and a band Dataset both feed a safe ast-based evaluator, which binds band tokens and evaluates each formula into an output stack of index layers. indices.yaml ndvi, evi, savi, ndwi … formula + bands band Dataset red, green, nir, swir … xarray DataArrays safe evaluator ast parse → validate bind tokens evaluate over arrays ndvi layer evi layer savi layer ndwi layer

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

Package Minimum version Why required
xarray 2023.6 DataArray container and vectorised arithmetic
rioxarray 0.15 Load GeoTIFF bands as georeferenced DataArrays
PyYAML 6.0 Parse the index registry
numpy 1.24 Backing array library for xarray 
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.

From a YAML entry to an evaluated array Each YAML block names an index, the bands it needs and the expression to evaluate. Validation checks that every referenced band exists in the loaded stack before anything runs, so a typo fails at parse time rather than three hours into a job. Only then is the expression evaluated against the chunked arrays. Fail at parse time, not at hour three indices.yaml - name: ndvi bands: [B08, B04] expr: (a-b)/(a+b) dtype: float32 validation every band present? expression parses? dtype and range sane? names unique? registry entry a callable bound to the named bands evaluated lazily, per chunk, once per index per tile Without the middle box a mistyped band name surfaces as a KeyError inside a worker, after the STAC search, after the reads, on tile 400 of 900 — with the partial outputs already written. 
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

One definition, three sensors Writing the index against semantic names such as red, nir and swir16, then supplying a per-sensor alias table, lets a single YAML entry run on Sentinel-2, Landsat 8 and MODIS. Without the alias layer each sensor needs its own copy of every definition, and they drift apart. Write the maths once, resolve the band names per sensor definition (nir − red) / (nir + red) alias table semantic name → this sensor's asset key Sentinel-2 L2A red → B04 · nir → B08 · swir16 → B11 Landsat 8/9 C2 L2 red → B4 · nir → B5 · swir16 → B6 MODIS MCD43A4 red → Band1 · nir → Band2 · swir16 → Band6 Band centres still differ between sensors — an alias makes the code run, not the values comparable. Record the sensor in the output metadata so a cross-sensor time series can be harmonised later.

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.

Frequently Asked Questions

Q: Why not just use Python eval for the index formulas? Python’s eval executes arbitrary code, so a malicious or malformed formula in the YAML could run anything. Parsing the formula with the ast module and walking only arithmetic nodes lets you evaluate math over DataArrays while rejecting function calls, attribute access, and imports.

Q: How does one engine handle indices with different bands? Each registry entry declares the band tokens it needs. The engine binds those tokens to DataArrays from the scene Dataset at evaluation time, so NDVI using red and nir and NDWI using green and nir both run through the same evaluator without special-casing.

Q: Can I add a new index without touching the code? Yes — add a new entry to the YAML registry with its name, formula, and bands. As long as the bands exist in the scene Dataset and the formula uses only arithmetic operators, the engine derives it with no code change.