Raster Subtraction with QgsRasterCalculator in PyQGIS

Question

I'm trying to do a simple subtraction of 2 raster layers using QgsRasterCalculator in PyQGIS. I've searched around and thought I nailed it but the code below didn't yield anything, but there wasn't even a crash.

What could have gone wrong?

Each raster only contains one band.

target_ras_input = 'D:/Users/abc/Desktop/target.tif'
target_ras = QgsRasterLayer(target_ras_input,"target_raster")
QgsMapLayerRegistry.instance().addMapLayer(target_ras)

initial_ras_input = 'D:/Users/abc/Desktop/initial.tif'
initial_ras=QgsRasterLayer(initial_bathy_input,"initial_raster")
QgsMapLayerRegistry.instance().addMapLayer(initial_ras)

target=QgsRasterCalculatorEntry()
target.raster = target_ras
target.bandNumber = 1
target.ref = 'target_raster@1'

initial=QgsRasterCalculatorEntry()
initial.raster = initial_ras
initial.bandNumber = 1
initial.ref = 'initial_raster@1'

entries = [ target , initial ]

final_output = 'D:/Users/abc/Desktop/final_raster.tif' 
calc=QgsRasterCalculator ( 'target.ref - initial.ref', final_output , 'GTiff', target_ras.extent(), target_ras.width(), target_ras.height(), entries )
calc.processCalculation()

final_raster=QgsRasterLayer(final_output, "please_work")
QgsMapLayerRegistry.instance().addMapLayer(final_raster)

Answer 1

Like Joseph said, there is an error in the first parameter of QgsRasterCalculator. But it must be string, unlike Josephs suggestion. So i tried this:

calc=QgsRasterCalculator ('%s - %s'%(target.ref, initial.ref), final_output , \
                      'GTiff', target_ras.extent(), target_ras.width(), \
                      target_ras.height(), entries )

You'll want to correct that line, too:

initial_ras=QgsRasterLayer(initial_bathy_input,"initial_raster")

because initial_bathy_input does nor exist, so change it to

initial_ras=QgsRasterLayer(initial_ras_input,"initial_raster")

Answer 2

I think you need to replace the single quotes in

'target.ref - initial.ref' 

otherwise the expression is read as a string and does not reference the actual names you used.

---

So it should read:

calc=QgsRasterCalculator( target.ref + ' - ' + initial.ref, final_output , 'GTiff', target_ras.extent(), target_ras.width(), target_ras.height(), entries )

Answer 3

For some who want the full working code with inputs as QgsRasterLayer:

def raster_substract(target_ras, initial_ras):
    target = QgsRasterCalculatorEntry()
    target.raster = target_ras
    target.bandNumber = 1
    target.ref = 'target_raster@1'
    initial = QgsRasterCalculatorEntry()
    initial.raster = initial_ras
    initial.bandNumber = 1
    initial.ref = 'initial_raster@1'
    entries = [target, initial]
    _, final_output = tempfile.mkstemp()
    formula = '%s - %s' % (target.ref, initial.ref)
    calc = QgsRasterCalculator(
        formulaString=formula,
        outputFile=final_output,
        outputFormat='GTiff',
        outputExtent=target_ras.extent(),
        nOutputColumns=target_ras.width(),
        nOutputRows=target_ras.height(),
        rasterEntries=entries)
    calc.processCalculation()
    return QgsRasterLayer(final_output, "DSM-DTM")