3

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

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Scott
  • 177
  • 1
  • 3
  • 12
  • 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. – Michael Stimson Jun 04 '14 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 Jun 04 '14 at 22:53
  • Here's a good link http://geoinformaticstutorial.blogspot.com.au/2012/09/reading-raster-data-with-python-and-gdal.html 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++. – Michael Stimson Jun 04 '14 at 23:02
  • Awesome, thank you! I will definitely look into that. – Scott Jun 04 '14 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 Jun 05 '14 at 00:24
  • @ChrisW, just the root of one.. there are other packages that are open source that could do the job that I'm not familiar with (or even aware of), if nothing better is presented it's at least a way out. – Michael Stimson Jun 05 '14 at 00:29

2 Answers2

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.

Community
  • 1
Chris W
  • 15,720
  • 2
  • 29
  • 47
1

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

#!/usr/bin/python 

# 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 
BandMath.ExecuteAndWriteOutput()

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

Community
  • 1
radouxju
  • 49,636
  • 2
  • 71
  • 144
  • 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 Jun 05 '14 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 Jun 05 '14 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 Jun 05 '14 at 23:14
  • User error of course, I forgot to add the ""s around the equation. It works, thanks again! – Scott Jun 05 '14 at 23:21