Adding all Binary Rasters in a Directory Using Raster Calculator in a Standalone Python Script

I see more code's been added since my previous answer :)

Ran this from the python console, you may need other imports if calling outside of QGIS.

Made a few small changes, namely

  • Added imports
  • need to set the QgsRasterCalculatorEntry.raster to rlayer (a QgsLayer instance), not the name of the layer
  • Added a list, layers, which gets filled with the layers in the loop
  • Added a one liner to create the complete sum expression.
  • Adds the result layer

This runs to completion. I'm afraid I don't have many good data sets to test this on. But it seems to work. As with your code it will use the extent and resolution of the first layer.

import os
import fnmatch
from PyQt4.QtCore import QFileInfo
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

inputFolder = "/path/to/rasters"
outputfilename = '/tmp/final.tif'

def findRasters(path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join(root, file)

entries = []
layers = []

for l in findRasters(inputFolder, '*.tif'):
    print l
    fileInfo = QFileInfo(l)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(l,baseName)
    QgsMapLayerRegistry.instance().addMapLayer(rlayer) 
    layer = QgsRasterCalculatorEntry()
    layer.ref = rlayer.name() + '@1'
    layer.raster = rlayer
    layers.append(rlayer)
    layer.bandNumber = 1
    print layer
    entries.append(layer) 

reflist = " + ".join([ent.ref for ent in entries]) 
expression = '(' +  reflist + ')'
print expression

calc = QgsRasterCalculator(expression,outputfilename,
                           'GTiff',
                           layers[0].extent(),
                           layers[0].width(),
                           layers[0].height(),
                           entries)

calc.processCalculation()
print "Finished!"

fileInfo = QFileInfo(outputfilename)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(outputfilename, baseName)
if not rlayer.isValid():
    print "Layer failed to load!"
else:
    QgsMapLayerRegistry.instance().addMapLayer(rlayer) 

You could use existing gdal tools.

there are two gdal command line utilities which you could call from your script using os.system or subprocess.call

gdalwarp - use this to create copies of each raster, scaled to the same size (use the -ts option to set resolution in pixels, or -tr in map units)

gdal_calc - this allows you to do the summation on the standardised rasters. If you have no more than 26 rasters you can do this in one call, otherwise you need to do this in a loop (add raster1+raster2, then add raster3, then add raster4 etc)