QGIS Script for Raster Equation

Question

I am trying to compute abs(A-B), where A is a multiple layer raster and B is a single raster.

I have tried the following script into the QGIS Processing Toolbox:

##XXX=name
##Utilities=group
##Raster=multiple raster
##B=raster
##OUT=folder

import glob, qgis
from PyQt4.QtCore import QFileInfo
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer

# Split rasters
layers = Raster.split(';')
output_path = OUT + "/"
suffix = "_suffix.tif"

# Get layer object
layer2 = processing.getObject(B)

for ras in layers:
    # Get layer object
    lyr = processing.getObjectFromUri(ras)
    entries = []    

    ras = QgsRasterCalculatorEntry()
    ras.ref = 'lyr@1'
    ras.raster = lyr
    ras.bandNumber = 1
    entries.append(ras)    

    raster2=QgsRasterCalculatorEntry()
    raster2.ref='B@1'
    raster2.raster=B
    raster2.bandNumber=1
    entries.append(raster2)

    calc = QgsRasterCalculator( 'abs(lyr@1-B@1)', output_path + lyr.name() + suffix, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
    calc.processCalculation()

for results in glob.glob(output_path + "*" + suffix): 
    fileInfo = QFileInfo(results)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

But it doesn't work. When I run the process I get the following message:

unicode cannot be converted to qgis._core.QgsRasterLayer in this context

Answer 1

1. Try replacing raster2.raster=B with raster2.raster=layer2.
2. Bizarrely, it seems you need to include spaces in your raster equation. So replace abs(lyr@1-B@1 with abs(lyr@1 - B@1).
3. I don't think you need to include the single raster into the same loop with the multiple rasters.

---

Here is the full code which ran successfully for me:

##XXX=name
##Utilities=group
##Raster=multiple raster
##B=raster
##OUT=folder

import glob, qgis
from PyQt4.QtCore import QFileInfo
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer

# Split rasters
layers = Raster.split(';')
output_path = OUT + "/"
suffix = "_suffix.tif"

# Get layer object
layer2 = processing.getObject(B)

for ras in layers:
    # Get layer object
    lyr = processing.getObjectFromUri(ras)
    entries = []        
    ras = QgsRasterCalculatorEntry()
    ras.ref = 'lyr@1'
    ras.raster = lyr
    ras.bandNumber = 1
    entries.append(ras)        

raster2=QgsRasterCalculatorEntry()
raster2.ref='B@1'    
raster2.raster=layer2
raster2.bandNumber=1
entries.append(raster2)

calc = QgsRasterCalculator( 'abs(lyr@1 - B@1)', output_path + lyr.name() + suffix, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
calc.processCalculation()

for results in glob.glob(output_path + "*" + suffix): 
    fileInfo = QFileInfo(results)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)