Evaluate spectral index expressions over raster bands and write a Cloud Optimized GeoTIFF — as a command line tool and as an importable Python library.
Every raster project grows its own ndvi.py. It hard-codes two band paths, loads both bands
fully into memory, divides them, and writes an untiled GeoTIFF. Then reality arrives: the red
band has nodata where the previous step clipped it, the scene has clouds that need the SCL band
masked out, the denominator hits zero over water and fills the output with inf, the DN values
need scaling to reflectance before EVI's coefficients mean anything, the file is 10 GB and no
longer fits in RAM, and the consumer downstream wants a COG with overviews.
band-math-cli is that script done properly, once. You give it a formula — either a built-in
index name or your own expression — and a mapping from band names to files. It validates that
the inputs actually share a grid, propagates every nodata source into the output mask, evaluates
the expression block by block, and writes a valid COG with statistics you can pipe into the next
stage.
The expression evaluator is deliberately not eval(). Expressions are parsed with Python's ast
module and every node is checked against an allowlist, so a formula coming from a config file, a
web form, or an untrusted job queue cannot reach the interpreter.
Python 3.12 or newer.
git clone https://github.com/python-remote-sensing/band-math-cli.git
cd band-math-cli
pip install -r requirements.txtThere is nothing to build or install — run it straight from the clone:
python -m band_math --helppython -m band_math --index ndvi \
--band red=B04.tif --band nir=B08.tif \
--scale 0.0001 --scl SCL.tif \
-o ndvi.tifwrote ndvi.tif
index : ndvi
expression : (nir - red) / (nir + red)
grid : 1098x1098 EPSG:32633
dtype : float32 (nodata nan)
overviews : 2, 4
stats : min=0.0534 max=0.8908 mean=0.5441 std=0.1779
valid : 986004/1205604 (81.8%)
--band name=file:N binds a name to band N of a multi-band file, so a single stack works
without splitting it first.
python -m band_math --index evi \
--band blue=stack.tif:1 --band red=stack.tif:3 --band nir=stack.tif:4 \
--scale 0.0001 --dtype int16 --json \
-o evi.tif{
"path": "evi.tif",
"expression": "2.5 * ((nir - red) / (nir + 6.0 * red - 7.5 * blue + 1.0))",
"index": "evi",
"bands": {
"blue": "stack.tif:1",
"red": "stack.tif:3",
"nir": "stack.tif:4"
},
"width": 1098,
"height": 1098,
"crs": "EPSG:32633",
"dtype": "int16",
"nodata": -32768.0,
"scale": 0.0001,
"offset": 0.0,
"compress": "deflate",
"blocksize": 512,
"overviews": [2, 4],
"stats": {
"valid_pixels": 1183644,
"total_pixels": 1205604,
"min": 0.0248823298686017,
"max": 1.18359624551575,
"mean": 0.4343757451094908,
"std": 0.2021766329069389,
"valid_fraction": 0.9817850637522769
}
}Any formula over your bound band names, using the allowlisted functions:
python -m band_math --expr "(nir - swir16) / (nir + swir16)" \
--band nir=B08.tif --band swir16=B11.tif \
--scale 0.0001 -o ndmi.tifUnsafe input is rejected before a single byte is read, with exit code 2:
python -m band_math --expr "__import__('os').system('id')" --band red=B04.tif -o x.tifband-math: only direct calls to allowlisted functions are permitted
python -m band_math --list-indices18 built-in indices:
arvi -- Atmospherically Resistant Vegetation Index
expression : (nir - (2.0 * red - blue)) / (nir + (2.0 * red - blue))
bands : blue, nir, red
range : [-1, 1]
reference : Kaufman & Tanre (1992), IEEE Trans. Geosci. Remote Sens. 30(2)
evi -- Enhanced Vegetation Index
expression : 2.5 * ((nir - red) / (nir + 6.0 * red - 7.5 * blue + 1.0))
bands : blue, nir, red
range : [-1, 1]
reference : Huete et al. (2002), Remote Sens. Environ. 83(1-2)
note : Requires surface reflectance in 0-1; DN input must be scaled first.
...
--list-indices --json emits the same registry as machine-readable JSON.
from band_math import INDEX_REGISTRY, compute_index
result = compute_index(
{"red": "B04.tif", "nir": "B08.tif"},
"ndvi.tif",
index="ndvi",
scale=0.0001,
scl="SCL.tif",
)
print(result.stats.mean, result.stats.valid_fraction)
print(INDEX_REGISTRY["ndvi"].reference)band_math.expr parses the expression with ast.parse(..., mode="eval") and walks every node
in the resulting tree. A node is rejected unless its type is on the allowlist: binary and unary
operators, comparisons, numeric literals, plain names, and direct calls to allowlisted functions.
That means Attribute, Subscript, Lambda, comprehensions, f-strings, dicts, lists, tuples
and string literals are all refused at parse time, which closes the whole family of sandbox
escapes that start with ().__class__.__bases__[0].__subclasses__(). Calls are additionally
checked: the callee must be a bare Name present in the function allowlist, keyword arguments
are refused, and arity is verified. Literal exponents above 64 are refused so red ** 1e9 can't
be used as a denial-of-service.
Evaluation is a small recursive walk that maps AST operator types to numpy ufuncs. There is no
eval, no compile, and no __builtins__ in scope at any point — the only things a name can
resolve to are the arrays you bound and the constants pi and e.
Available functions: abs, clip, exp, log, log10, maximum, minimum, power,
sign, sqrt, square, where.
Rasters are opened once per distinct path and compared against the first one: CRS must be equal,
pixel dimensions must be equal, and the six affine coefficients must match within 1e-6. If they
don't, you get a diagnostic naming the offending file and the specific difference, rather than a
silently misaligned result:
band-math: input rasters do not share a common grid (relative to B04.tif):
B11_20m.tif: shape 549x549 != 1098x1098
Reproject or resample the inputs onto a single grid first.
An --scl or --mask raster is checked against the same grid.
Four things can invalidate a pixel, and all four are unioned into one mask:
- Per-band nodata. Each band is read with
read_masks(), so both the declared nodata value and any internal per-dataset mask band are honoured. A pixel invalid in any input band is invalid in the output. - Non-finite results. Division by zero (
nir + red == 0over deep water or a burn scar) and out-of-domain math (sqrtof MSAVI2's negative discriminant) are evaluated undernp.errstate(all="ignore"), then the result is tested withnp.isfinite. Anything non-finite becomes nodata. You never getinfor a straynanwritten as if it were data, and numpy never emits a warning storm mid-pipeline. - SCL classes.
--scl SCL.tifdrops classes0,1,2,3,8,9,10,11by default — no-data, defective, cast shadows, cloud shadows, medium and high probability cloud, thin cirrus, and snow/ice. Vegetation, bare soil, water and unclassified are kept. Override with--scl-drop, using numbers or names (--scl-drop cloud_shadows,cloud_high_probability). - An arbitrary mask raster.
--mask keep.tifkeeps pixels where the mask is non-zero.
Optionally, --drop-out-of-range also masks results outside a built-in index's documented valid
range, and --clip-range clamps them into it instead.
Sentinel-2 L2A distributes reflectance as scaled integers. --scale 0.0001 converts DN to
reflectance before evaluation; from processing baseline 04.00 onwards you also need
--offset -0.1. This matters more than it looks: NDVI is a ratio and survives an unapplied
scale, but EVI's 6.0 * red - 7.5 * blue + 1.0 has an absolute +1 in it, so feeding it raw DN
gives numerically meaningless output that still looks like a plausible raster.
On the output side, --dtype int16 stores the index multiplied by --out-scale (default 10000),
which roughly halves file size versus float32 while keeping four decimal places. The scale factor
is written into the GeoTIFF so readers can recover physical values.
The output grid is tiled into windows (--chunk, default max(blocksize, 512)) and each window
is read, evaluated, masked and written independently, so peak memory is a function of the chunk
size and band count, not of the file size. Statistics are accumulated in streaming fashion
(count, sum, sum of squares, running min/max), so the JSON summary is exact regardless of chunk
size — the test suite asserts that changing the chunk size changes neither a single output pixel
nor a single statistic.
Blocks are written to a temporary tiled GeoTIFF next to the destination, then handed to GDAL's
COG driver via rasterio.shutil.copy, which lays out the IFDs in the required order and builds
internal overviews down to the block size. The result passes COG validation: tiled, overviews
present, headers at the front of the file.
| Flag | Default | Description |
|---|---|---|
--index NAME |
— | Built-in index to compute (see --list-indices) |
--expr EXPRESSION |
— | Custom expression over the bound band names |
--list-indices |
— | Print the registry (add --json for JSON) and exit |
--band NAME=PATH[:N] |
— | Bind a band name to a file, or to band N of it. Repeatable |
-o, --output PATH |
— | Output COG path |
-f, --overwrite |
off | Replace the output if it exists |
--scl PATH |
— | Sentinel-2 Scene Classification band to mask with |
--scl-drop CLASSES |
0,1,2,3,8,9,10,11 |
SCL classes to drop, as numbers or names |
--mask PATH |
— | Boolean mask raster; non-zero pixels are kept |
--clip-range |
off | Clamp results into the index's valid range |
--drop-out-of-range |
off | Mask results outside the index's valid range |
--scale F |
1.0 |
Multiply every input band by F before evaluation |
--offset F |
0.0 |
Add F to every input band after scaling |
--dtype T |
float32 |
float32, float64, int16 or int32 |
--nodata V |
per dtype | Output nodata (nan for float, -32768 for int16) |
--out-scale F |
10000 |
Storage scale for integer dtypes |
--compress C |
deflate |
deflate, zstd, lzw or none |
--blocksize N |
512 |
COG internal tile size (multiple of 16) |
--chunk N |
max(blocksize, 512) |
Processing window size in pixels |
--json |
off | Print a JSON summary on stdout |
-q, --quiet |
off | Suppress the text summary |
Exit codes: 0 success, 1 computation failure, 2 usage or validation error.
| Object | Purpose |
|---|---|
compute_index(bands, output, *, index=..., expression=..., ...) |
Full pipeline; returns an IndexResult |
compute_array(expression, bands, valid=None, *, scale=1.0, offset=0.0) |
Evaluate over in-memory arrays; returns (values, valid) |
safe_eval(expression, variables) |
Parse, validate and evaluate an expression |
parse_expression(expression) |
Validate only; raises ExpressionError |
expression_variables(expression) |
The band names an expression requires |
INDEX_REGISTRY / get_index(name) |
The built-in SpectralIndex definitions |
IndexResult / Stats |
Result metadata; .to_dict() for JSON |
BandSpec / parse_band_spec / BandStack |
Band binding and grid validation |
SCL_CLASSES / parse_scl_classes / scl_valid_mask |
SCL masking helpers |
Exceptions: ExpressionError, BandBindingError, GridMismatchError, MaskError,
BandMathError.
Built-in indices use sensor-neutral names — blue, green, red, rededge, nir, swir16,
swir22 — which map to Sentinel-2 B02/B03/B04/B05/B08/B11/B12 and Landsat 8/9
B2/B3/B4/-/B5/B6/B7. Custom expressions can use any identifier you like.
- No reprojection or resampling. Inputs must already be on a common grid; the tool tells you precisely what mismatches rather than guessing. Resample Sentinel-2 20 m bands to 10 m first if you want to mix them.
- Single-band output. One expression, one band, one file. Computing several indices means several invocations (cheap — they share the OS page cache for the inputs).
- No temporal or multi-scene logic. No compositing, no mosaicking, no time series.
- Overviews only appear when the raster is larger than the block size. A 256×256 output with the default 512 block size is a valid COG with no overview levels.
- Statistics are computed on the physical index values, before integer storage scaling, and reflect only unmasked pixels.
--scale/--offsetapply to every band uniformly. Mixing bands with different scale factors requires pre-scaling them yourself.- Grid comparison is exact (1e-6 on the affine coefficients). Rasters that are conceptually aligned but carry floating-point drift in their transforms will be rejected.
Background on the concepts this tool implements, from the guides at Python Remote Sensing:
- Computing EVI and NDWI from Sentinel-2 bands walks through the coefficients and the reflectance scaling that these formulas assume.
- Masking clouds with the Sentinel-2 SCL band covers what each SCL class actually means and which ones are worth dropping.
- Understanding Cloud Optimized GeoTIFF structure explains the IFD ordering and overview layout that make the output here range-readable.
- Rasterio vs xarray for band math is a useful comparison if you are deciding whether a windowed rasterio tool like this one or a labelled-array approach fits your pipeline better.
- Optimizing rasterio window reads for memory efficiency goes deeper on choosing chunk sizes for out-of-core processing.
Issues and pull requests are welcome. To work on the code:
pip install -r requirements-dev.txt
ruff check . && ruff format --check .
python -m pytestNew built-in indices need an entry in band_math/indices.py with a literature reference and a
valid range; the test suite checks every registered index parses, declares the bands its formula
uses, and produces an in-range value for plausible reflectance input.
MIT — see LICENSE.
Built and maintained alongside Python Remote Sensing & Raster Processing Pipelines.