gdal_calc

What is gdal_calc?

gdal_calc.py is GDAL's command-line raster calculator. It binds input rasters to single-letter variables (A, B, C, …), evaluates a Python/NumPy expression pixel-by-pixel, and writes the result to an output raster. It is how you compute indices like NDVI, apply thresholds, reclassify, or combine rasters arithmetically without writing a Python script.

gdal_calc.py -A <input1> [-B <input2>...] --outfile=<out.tif> --calc="<expression>"

Commonly used options:

• -A, -B, … — input rasters (any OGR-opened letter)
• --A_band <n> — source band for each letter (default 1)
• --calc="<expression>" — NumPy expression (supports operators, logical_and, where, minimum, etc.)
• --outfile=<filename> — output raster
• --NoDataValue=<value> — output NoData
• --type=<Byte|Int16|Float32...> — output data type
• --format=<driver> / --co=<NAME>=<VALUE> — output driver and creation options
• --allBands=A — evaluate the expression per band of A, producing a multi-band output
• --overwrite — overwrite existing output
• --quiet / --debug — verbosity

When would you use gdal_calc?

Use gdal_calc.py for any pixel-wise arithmetic that does not need a neighbourhood operation. Typical jobs: computing NDVI (gdal_calc.py -A nir.tif -B red.tif --calc="(A.astype(float)-B)/(A+B)" --outfile=ndvi.tif --type=Float32 --NoDataValue=-9999), thresholding an index to a binary mask (--calc="A>0.3"), reclassifying values (--calc="where(A==1, 10, where(A==2, 20, 0))"), or summing two rasters with matching grids.

For operations that need a moving window (slope, focal statistics), use gdaldem or a dedicated tool. For very large rasters, gdal_calc.py is convenient but not the fastest — it reads bands into memory block-wise. When performance matters, consider writing a short Python script using GDAL's ReadAsArray with explicit block iteration, or use numexpr-backed tools like rio calc.

FAQs

Why is my NDVI output all zeros or all ones?

Integer division. If A and B are UInt16, (A-B)/(A+B) is integer-divided and quantised to 0 or 1. Cast to float inside the expression: (A.astype(float)-B)/(A+B). Also set --type=Float32 on the output; otherwise the result is truncated to the default Float32-or-source type which may also surprise you. Finally, verify NoData handling: Sentinel-2 saturates to 0 outside the swath, which will divide by zero.

How do I reclassify with gdal_calc?

Use numpy.where nested calls: --calc="where(A<100, 1, where(A<200, 2, where(A<500, 3, 4)))". For many classes, this gets clumsy; gdal_translate with a VRT colour map, gdalwarp -r mode, or a dedicated reclass tool like whitebox_tools or rio are more legible. For table-driven reclass, pre-write a Python script.

Can I evaluate expressions across all bands in one input?

Yes — use --allBands=A. The expression is evaluated per band, producing an output with the same band count as A. Useful for applying the same scaling or mask to every band of a multispectral scene in one pass.

How do I handle NoData correctly?

Pass --NoDataValue to set the output NoData, and make sure the input NoData values are also declared on the source rasters (check with gdalinfo). gdal_calc.py masks NoData input pixels so they do not participate in arithmetic, but only if the source advertises a NoData value. If NoData is implied via a sentinel like -9999 but not set, use gdal_edit -a_nodata -9999 to declare it first.