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.– Michael StimsonJun 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.– ScottJun 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++.– Michael StimsonJun 4, 2014 at 23:02
-
Awesome, thank you! I will definitely look into that.– ScottJun 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 WJun 5, 2014 at 0:24
|
Show 1 more comment
2 Answers
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:
- What is the equivalent of arcpy "Con" in QGIS and/or R raster-package?
- Replacing raster values in one raster with another raster in QGIS ... how to?
- QGIS Raster Calculater Syntax
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.
#!/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
-
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!– ScottJun 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?– ScottJun 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 > ?– radouxjuJun 5, 2014 at 23:14
-
User error of course, I forgot to add the ""s around the equation. It works, thanks again!– ScottJun 5, 2014 at 23:21