I'm wondering if there is a method to perform a function similar to Arc's Con tool, using open source options in python. I have two rasters, one that shows locations of clouds and the other is just a landsat band. If a value is >x in the cloud raster, I want to take the pixel from the landsat band, and if it's

  • It's something you could do easily with GDAL, but you'd have to do the hard work yourself, reading each image to a buffer, comparing buffers and writing an output buffer to a new image. Commented Jun 4, 2014 at 22:41
  • Would you mind going into a little more detail? I'm not a complete python noob, but I wouldn't say I'm an expert either. For the research project I'm working on, open source is fairly important so I'd like to try and figure out how to do this, rather than resort to the arcpy module.
    – Scott
    Commented Jun 4, 2014 at 22:53
  • Here's a good link geoinformaticstutorial.blogspot.com.au/2012/09/… if the cells line up there's no problem, the fun starts if the rasters are different size or origin. Generally I don't use GDAL in python, but have done; most of my work with GDAL is in C++. Commented Jun 4, 2014 at 23:02
  • Awesome, thank you! I will definitely look into that.
    – Scott
    Commented Jun 4, 2014 at 23:30
  • Your answer appears cut off, so some specific approaches might be left off. @MichaelMiles-Stimson I'd say you've got an answer there or at least the root of one.
    – Chris W
    Commented Jun 5, 2014 at 0:24

2 Answers 2


Con is a conditional raster calculation, therefore it could be done as a series of steps or directly with a raster calculator/functions that support conditionals. There are several options already covered at other questions here on SE. This may actually be a duplicate of one, but I will just link them since you are asking for any python based open source option and these are specific to QGIS, SAGA, and R:

As mentioned in some of those answers, QGIS with the RasterCalc plugin (not the default raster calculator) supports conditional statements.


apart from QGIS, you can also use OTB, which implements the muParser functions and can be wrapped in Python. It is very efficient.


# Import the otb applications package 
import otbApplication 

# The following line creates an instance of the BandMath application 
BandMath = otbApplication.Registry.CreateApplication("BandMath") 

# The following lines set all the application parameters: 
BandMath.SetParameterStringList("il", ['cloud_r.tif', 'landsat_nir.tif']) 

BandMath.SetParameterString("out", "apTvUtBandMathOutput.tif") 

BandMath.SetParameterString("exp", "(im1b1>1)?(im2b1):(0)") 

# The following line execute the application 

Another alternative is gdal_calc.py, which can be used like here (can be launched with subprocess.call() from Python )

gdal_calc.py -A cloud.tif -B landsat_NIR.tif --outfile=out_NIR.tif --calc="B*(A>1)" 

note: I assume here that your "x" is 1 , but you can change the value

  • I like this option, I've used the gdal_calc.py file once before so I'm already familiar with it. I had no idea I could use it like the con tool though. Thanks!
    – Scott
    Commented Jun 5, 2014 at 19:48
  • I tried to use this method, but the > is causing problems. Does this not work if I'm calling the script inside of another python script with the os.system command? I.e. can you only call this from command line?
    – Scott
    Commented Jun 5, 2014 at 23:04
  • it should work but you'll need to call it with "python gdal_calc.py ...". Note that I prefer subprocess.call instead of os.system. What's wrong with the > ?
    – radouxju
    Commented Jun 5, 2014 at 23:14
  • User error of course, I forgot to add the ""s around the equation. It works, thanks again!
    – Scott
    Commented Jun 5, 2014 at 23:21

Your Answer

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

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