Combining two rasters based on a Condition using Open Source Python

Question

I'm new to Python. My only experience is with Arcpy. I need to combine two rasters with the same extent into one output raster based on a condition. Using Arcpy, I would do something like this:

outCon = Con ((inRaster1 > 0.5), inRaster2, inRaster1)

but I don't have access to Arcpy and need to use an Open Source Python (Rasterio, numpy, gdal, ...). Any ideas or thoughts?

Answer 1

GDAL_Calc should be able to do what you're after, but from python it would need to be shelled with os.system or subprocess.Popen as it is a python script itself:

GDAL_Calc.py -A inRaster1 -B inRaster2 --calc="((A > 0.5) * B) + ((A <= 0.5) * A)" ... then the rest of the args depending on what your intention is

The condition (A > 0.5) is 1 if inRaster1 is greater than 0.5 making the condition (A <= 0.5) equaling 0, making the sum ((1)*B)+((0)*A), conversely if A <= 0.5 the conditions switch making the sum ((0)*B)+((1)*A).

Be sure to run from an OSGeo or GDAL shell (two options of obtaining GDAL on Windows) so that the script can find the libs.

Answer 2

Thank you for the answer. I decided to give np.where a go. I tried the following:

import os
import rasterio
import numpy as np

workspace = '.../tests'
A = os.path.join (workspace, "inRaster1.tif")
B = os.path.join (workspace, "inRaster2.tif")

with rasterio.open(A) as srcA:
     array_inA = srcA.read()
     profile = srcA.profile

     with rasterio.open(B) as srcB:     
        array_inB = srcB.read()

        array_out= np.copy(array_inA)
        array_out[np.where(array_inA > 0.5),array_inB ,array_inA]

result = os.path.join (workspace, "result.tif")  
with rasterio.open(result, 'w', **profile) as dst:
    dst.write(array_out)
    crs = srcA.crs

But I'm getting this error.

line 29, in <module> array_out[np.where(array_inA< 0.5),array_inB ,array_inA]
IndexError: arrays used as indices must be of integer (or boolean) type