21

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?

1

4 Answers 4

17

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

1
  • 1
    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 5, 2016 at 19:47
14

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
1

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
1

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'
     })

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.