0

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?

2

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.

|improve this answer|||||
0

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
|improve this answer|||||
  • This would probably be better suited as an edit to the original post since it's not an answer per se. But it looks like an error in syntax, the line to create the output should be array_out = np.where((array_inA > 0.5), array_inB, array_inA). Note that the 2 arrays are being provided as arguments to np.where. If they're not, where returns the indexes at which the condition is satisfied. And there's no need to make a copy first. – mikewatt Feb 20 at 19:18
  • The code works now. Thanks. Sorry, it is my first post here. I will make sure to get familiar with the stack exchange rules. – MoeGIS Feb 20 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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