Batch Computing Indices Across a STAC Collection

To run an index pipeline over an entire search, freeze the item list, make each item independently processable, and record the outcome of each:

items = list(search.items())
manifest = [{"id": i.id, "datetime": i.datetime.isoformat(), "tile": i.properties.get("s2:mgrs_tile")}
            for i in items]

results = [process_item(entry) for entry in manifest]      # or a pool over the same list

This is the operational form of the pipeline described in Spectral Index Calculation Pipelines.

Why This Arises in Remote Sensing Workflows

An index pipeline that works on one scene rarely works on ten thousand without change. The arithmetic is unchanged; what breaks is everything around it. One scene is missing an asset. One has a null CRS. One is 403 because a token expired mid-run. A worker is preempted. The process dies at 80 percent and nobody knows which scenes completed.

The structural answer is to make the item the unit of work and to make that unit idempotent: a function of an item identifier that produces a deterministic output and can be run twice with no harm. Every operational property follows from that — retry, resume, parallelism, partial delivery — and none of them are available if the pipeline is written as a loop that accumulates state.

One item in, one of three outcomes out Each manifest entry is processed independently. If its output already exists and validates it is skipped; if processing succeeds a new output and a status row are written; if it fails the error is recorded and the item is quarantined. No outcome affects any other item. Independent units, recorded outcomes manifest row item id + date output exists? deterministic path + validates? yes no read → mask → index → write atomic: temp then rename written status: ok failed error recorded skipped, already done Re-running the whole manifest is safe — that is what makes resume trivial.

Environment & Setup

Package Version Why
pystac-client ≥0.7 Building the work list
rioxarray ≥0.15 Reads, index arithmetic, COG writes
rasterio ≥1.3.0 Validation of existing outputs
pandas ≥2.0 Manifest and status table 
pip install "pystac-client>=0.7" "rioxarray>=0.15" "rasterio>=1.3.0" "pandas>=2.0"

Complete Working Example

This module builds a frozen manifest, processes one item idempotently, and runs the whole list with failures isolated.

Per-item statuses a batch run should distinguish A batch run has more outcomes than success and failure. Distinguishing skipped, written, transient failure, permanent failure and no-input tells you which items to retry, which to investigate and which are correctly absent. Five outcomes, not two status meaning what to do written output produced this run nothing skipped output already existed and validated nothing — this is resume working transient failure timeout, 503, throttling re-run the same manifest permanent failure missing asset, bad CRS investigate; retrying will not help no input search returned no assets expected for some tiles and dates Lumping the last three together is what makes a failed batch run hard to recover from. 
import os
import tempfile
from concurrent.futures import ThreadPoolExecutor, as_completed

import numpy as np
import pandas as pd
import rasterio
import rioxarray  # noqa: F401
import xarray as xr
from rasterio.enums import Resampling

INDEX_DEFS = {
    "ndvi": ("nir", "red"),
    "ndwi": ("green", "nir"),
    "nbr": ("nir", "swir22"),
}
ASSETS = {"red": "B04", "green": "B03", "nir": "B08", "swir22": "B12", "scl": "SCL"}


def build_manifest(search, path: str) -> pd.DataFrame:
    """Freeze the work list so the run is reproducible and resumable."""
    rows = [{"id": i.id,
             "datetime": i.datetime.isoformat(),
             "tile": i.properties.get("s2:mgrs_tile"),
             "cloud": i.properties.get("eo:cloud_cover")}
            for i in search.items()]
    df = pd.DataFrame(rows).sort_values("datetime")
    df.to_parquet(path, index=False)
    return df


def output_path(root: str, item_id: str, index: str) -> str:
    """Deterministic: the same item always maps to the same file."""
    return os.path.join(root, index, f"{item_id}_{index}.tif")


def already_done(path: str) -> bool:
    """Skip only if the existing output is readable and complete."""
    if not os.path.exists(path):
        return False
    try:
        with rasterio.open(path) as src:
            return src.width > 0 and src.count >= 1
    except Exception:
        return False                       # corrupt leftovers are reprocessed


def nd(a: xr.DataArray, b: xr.DataArray) -> xr.DataArray:
    denom = a + b
    return xr.where(denom != 0, (a - b) / denom, np.nan).astype("float32")


def process_item(item, root: str, indices=("ndvi",)) -> dict:
    """Compute the requested indices for one item. Never raises; returns a status row."""
    status = {"id": item.id, "written": [], "skipped": [], "error": None}
    try:
        wanted = [ix for ix in indices
                  if not already_done(output_path(root, item.id, ix))]
        status["skipped"] = [ix for ix in indices if ix not in wanted]
        if not wanted:
            return status

        needed = {b for ix in wanted for b in INDEX_DEFS[ix]} | {"scl"}
        bands = {}
        for key in needed:
            href = item.assets[ASSETS[key]].href
            da = rioxarray.open_rasterio(href, masked=True, chunks={"x": 1024, "y": 1024})
            bands[key] = da.squeeze(drop=True).astype("float32")

        reference = bands[INDEX_DEFS[wanted[0]][0]]
        for key, da in bands.items():
            if da.rio.transform() != reference.rio.transform():
                method = Resampling.nearest if key == "scl" else Resampling.bilinear
                bands[key] = da.rio.reproject_match(reference, resampling=method)

        clear = ~bands["scl"].isin([0, 1, 3, 8, 9, 10, 11])

        for ix in wanted:
            a_key, b_key = INDEX_DEFS[ix]
            arr = nd(bands[a_key].where(clear), bands[b_key].where(clear))
            arr.attrs.update(index=ix, item_id=item.id, datetime=item.datetime.isoformat())

            dst = output_path(root, item.id, ix)
            os.makedirs(os.path.dirname(dst), exist_ok=True)
            fd, tmp = tempfile.mkstemp(suffix=".tif", dir=os.path.dirname(dst))
            os.close(fd)
            try:
                arr.rio.to_raster(tmp, driver="COG", compress="DEFLATE", blocksize=512)
                os.replace(tmp, dst)          # atomic: a killed worker leaves no partial output
            finally:
                if os.path.exists(tmp):
                    os.remove(tmp)
            status["written"].append(ix)
    except Exception as exc:
        status["error"] = repr(exc)
    return status


def run(items, root: str, indices=("ndvi", "ndwi"), workers: int = 8) -> pd.DataFrame:
    """Process every item, isolating failures."""
    rows = []
    with ThreadPoolExecutor(max_workers=workers) as pool:
        futures = [pool.submit(process_item, it, root, indices) for it in items]
        for fut in as_completed(futures):
            rows.append(fut.result())
    df = pd.DataFrame(rows)
    df.to_parquet(os.path.join(root, "status.parquet"), index=False)
    return df


if __name__ == "__main__":
    status = run(list(search.items()), root="out/indices")
    print("ok:", int(status["error"].isna().sum()), "failed:", int(status["error"].notna().sum()))
    print(status[status["error"].notna()].head())

Three properties make this operationally sound. The output path is a pure function of the item and the index, so nothing depends on run order. The write is atomic through a temporary file and os.replace, so a killed worker never leaves a half-written COG that the next run mistakes for a completed one. And process_item never raises, so one broken scene costs one row in the status table rather than the whole job.

Variant Patterns

1. Threads, processes or a cluster

Pick the executor from the bottleneck An index batch over remote COGs spends most of its time waiting on the network, so threads are the right executor. Once the arithmetic is heavier — many indices per scene, or resampling — processes win. Only when a single unit exceeds one machine's memory does a distributed cluster earn its complexity. Executor by bottleneck, not by preference threads reads dominate GDAL releases the GIL 8–16 workers the default for this job processes arithmetic dominates many indices per scene workers ≈ cores watch memory per worker cluster one unit exceeds a machine or a cross-scene reduction chunked arrays, lazy graph complexity worth paying for A per-scene index batch belongs in the left column far more often than teams assume.

For the distributed case the same process_item becomes a task submitted to a scheduler, as in Computing NDVI with xarray on a Dask Cluster — the function does not change, only who calls it.

2. Externalising the index definitions

Hard-coding INDEX_DEFS works until someone wants a new index. Moving them into configuration turns that into an edit rather than a deployment, which is the argument made in Building a YAML-Driven Multi-Index Pipeline.

3. Streaming a summary alongside the rasters

Most consumers of a batch index run want statistics, not imagery. Computing them during the same pass costs almost nothing extra, because the array is already in memory.

summary = {
    "id": item.id, "index": ix,
    "valid_fraction": float(arr.notnull().mean()),
    "median": float(arr.median(skipna=True)),
    "p10": float(arr.quantile(0.1, skipna=True)),
    "p90": float(arr.quantile(0.9, skipna=True)),
}

Written to a table, these become a time series that can be plotted without touching a single raster again — the same shape recommended in Zonal Statistics and Vector–Raster Integration.

Monitoring a Long Run

Three numbers tell you whether a batch is healthy, and all three come from the status table.

The failure rate should be low and its causes few. A handful of failures with the same exception is a bug worth fixing; hundreds of distinct exceptions usually means credentials or network trouble rather than data trouble.

The skip rate tells you whether the run is doing what you think. On a resumed run it should start high and fall; on a first run it should be zero, and a non-zero value means outputs from an earlier configuration are being silently accepted.

Throughput, measured as items per minute, should be roughly flat. A steady decline usually means memory pressure building in a long-lived process, and a sudden collapse means throttling — the bucket asking for fewer requests, which is handled by backing off rather than by retrying harder, as covered in Retrying Failed Raster Tasks in a Prefect Pipeline.

Common Errors

The run stops on a single bad scene

An exception escaped the per-item function. Catch broadly inside process_item and record the error; only genuinely fatal conditions — bad credentials, missing output directory — should stop the run.

A resumed run reprocesses everything

The output path depends on something that changes between runs, such as a timestamp. Derive it from the item id alone.

Outputs exist but are unreadable

Writes were not atomic and a worker was killed mid-write. Write to a temporary file and os.replace it into place.

Frequently Asked Questions

Q: Should the unit of work be a scene or a tile? A scene, when indices are computed per acquisition and written per acquisition. It matches the STAC item, needs no coordination, and makes retry trivial. Split into sub-tiles only when a single scene exceeds a worker’s memory.

Q: Why write a manifest instead of re-running the search? Because a catalogue changes. Items get added, reprocessed or withdrawn, so two searches days apart return different sets and the run stops being reproducible. Freeze the list, then process it.

Q: How do I resume a run that died halfway? If each output path is deterministic and each write is atomic, resuming is simply re-running: completed items are skipped because their outputs exist and validate.