gdal_calc raster calculator syntax for logical operators and other functions

Question

In documentation for gdal_calc it is stated Command line raster calculator with numpy syntax. Later on there are few examples where in one of them:

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0 - means set values of zero and below to null

Unfortunately there is no example on logical operators like:

--calc="A*(A>0 and A>B)"- means keep A if A bigger zero and bigger B and set the rest to null

Based on Numpy/Scipy logic functions I would expect to write logical operators as:

--calc="A*logical_and(A>0,A>B)"

I tried this and it seems to work but I would like to be assured that is correct.

In the similar way if you want minimum of A and B:

--calc="A*(A<=B)+B*(A>B)"

You can just write:

--calc="minimum(A,B)"

My issue is I can't find any cookbook to make sure I get this right. Is there some good cookbook with advanced examples of what is and is not possible with gdal_calc?

--calc=where(condition, if, else) ? See https://gis.stackexchange.com/q/238397/6998 — Campa, Dec 06 '22 at 18:23

score 17 · Accepted Answer · edited Dec 04 '17 at 12:10

In the source for gdal_calc.py, the calculation is made directly using eval:

myResult = eval(opts.calc, global_namespace, local_namespace)

That would suggest that any well-formed expression that also evaluates on the command line will work. According to the documentation, you may use gdalnumeric syntax with +-/*, and/or numpy functions. You can test your functions using small dummy arrays in the interactive shell, then use the same calls in gdal_calc.

Keep in mind that chaining together multiple numpy functions is likely to produce temporary in-memory arrays that can substantially increase memory usage, especially when dealing with large images.

You can look at the numpy documentation for a list of all the functions: routines. The ones you are after are likely here: math or here: routines.logic.

This is where functions like minimum are coming from, it's just that the namespace is already imported. Really, it's numpy.minimum, etc

Thank you Ben, that is another way I did not have a clue about. Still after some cookbook which would explain what is possible to use, because eval does not include minimum() etc functions which are actually possible to use in expression. — Miro, Aug 05 '16 at 19:47

score 14 · Answer 2 · answered May 27 '17 at 02:01

Following on from Benjamin's answer, you can use logical_or() or logical_and(). See http://docs.scipy.org/doc/numpy/reference/routines.logic.html. The following example worked nicely for me. This sets all values between 177 and 185 (inclusive) to 0, which is then treated as nodata.

gdal_calc.py -A input.tif --outfile=output.tif --calc="A*logical_or(A<=177,A>=185)" --NoDataValue=0

score 1 · Answer 3 · answered Feb 16 '18 at 15:58

I had a raster where values ranged between -1 and 3 where zero is a valid number. I had some problems making a gdal_calc expression so made this fast and furious solution.

#!/usr/bin/env python3

fileNameIn = "/tmp/geotiff/Global_taxonomic_richness_of_soil_fungi.tif"
fileNameOut = "/tmp/geotiff/Global_taxonomic_richness_of_soil_fungi.tiff"
dst_options = ['COMPRESS=DEFLATE',"PREDICTOR=3","TILED=YES"]
noDataValue = -3.4028234663852886e+38

from osgeo import gdal
import numpy

src_ds = gdal.Open(fileNameIn)
format = "GTiff"
driver = gdal.GetDriverByName(format)
dst_ds = driver.CreateCopy(fileNameOut, src_ds, False ,dst_options)

# Set location
dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
# Set projection
dst_ds.SetProjection(src_ds.GetProjection())
srcband = src_ds.GetRasterBand(1)

dataraster = srcband.ReadAsArray().astype(numpy.float)
#Rplace the nan value with the predefiend noDataValue
dataraster[numpy.isnan(dataraster)]=noDataValue

dst_ds.GetRasterBand(1).WriteArray(dataraster)
dst_ds.GetRasterBand(1).SetNoDataValue(noDataValue)

dst_ds = None

score 1 · Answer 4 · edited Sep 08 '21 at 11:35

This is yet another way that works to do it. I'm using GDAL via QGIS in PyQGIS. Anything can be specified in the 'FORMULA'.

surfaceA = 'P:/70210 Severn Trent Groundwater Map/GIS/OS_100k_WaterFeatures_Tamworth_DTM_interpolated_50mRes_WOCanal.tif'
surfaceB = 'P:/70210 Severn Trent Groundwater Map/GIS/Main_river_points_DTM_inTamworth_interpolated_50mRes.tif',
processing.run("gdal:rastercalculator", 
    {'INPUT_A': surfaceA,
     'BAND_A':1,
     'INPUT_B': surfaceB, 
     'BAND_B':1,
     'FORMULA':'numpy.min((A,B),axis=0)',
     'NO_DATA':None,
     'RTYPE':5,
     'OPTIONS':'',
     'EXTRA':'',
     'OUTPUT':'TEMPORARY_OUTPUT'
     })

gdal_calc raster calculator syntax for logical operators and other functions

4 Answers4

Linked