2

I'm writing a QGIS plugin and I need to create a raster layer as the average of other rasters. How can I do it? I have tried to use gdal_calc.py through os.system, something like:

cmd = 'gdal_calc.py 
       -A first_layer.tif 
       -B second_layer.tif 
       --outfile=AVG.tif 
       --cal="(A+B)/2"'
os.system(cmd)

but it seems that the command (which works perfectly fine in the command line) is not executed from my QGIS plugin. What should I use to create my average raster?


I finally found a solution, and I share it for those who might be interested:

from gdal_calculations import Dataset, Env


def aggregate_layers(self, layers):
    datasets = []
    for l in layers:
        datasets.append(Dataset(l))
    sum = datasets[0]
    for i in range(1, len(datasets) - 1):
        sum += datasets[i]
    Env.overwrite = True
    avg = sum
    avg /= len(datasets)
    root = '/home/kalimaha'
    file_name = root + '/average.tif'
    avg.save(file_name)
    if self.dlg.open_in_qgis.isChecked() is True:
        rl = self.iface.addRasterLayer(file_name, 'Average Layer')

I have tried @xunilk's code, but my QGIS simply crashes. I had to edit it a bit, this is the final code:

def aggregate_layers(self, layers):

    entries = []
    for l in layers:
        fileInfo = QFileInfo(l)
        baseName = fileInfo.baseName()
        rlayer = QgsRasterLayer(l, baseName)
        QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')
        layer = QgsRasterCalculatorEntry()
        layer.ref = rlayer.name() +'@1'
        layer.raster = rlayer
        layer.bandNumber = 1
        entries.append(layer)

    for e in entries:
        QgsMessageLog.logMessage('e: ' + str(e), 'Geobricks TRMM')

    fileInfo = QFileInfo(layers[0])
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(layers[0], baseName)
    QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')

    expression = '(' + entries[1].ref + ' + ' + entries[0].ref + ') / 2'
    QgsMessageLog.logMessage('Expression: ' + str(expression), 'Geobricks TRMM')
    calc = QgsRasterCalculator(expression,
                               '/Development/GIS/average.tif',
                               'GTiff',
                               rlayer.extent(),
                               rlayer.width(),
                               rlayer.height(),
                               entries)
    QgsMessageLog.logMessage('calc: ' + str(calc), 'Geobricks TRMM')

My layers are not in QGIS yet, so I have to load it from the filesystem, the array nemd layers contains the absolute paths of the layers I need to aggregate. So, when I do calc = QgsRasterCalculator(... QGIS crashes.

6
  • Try putting python or python.exe in front of your gdal_calc.py or use a full path. Those might help.
    – Branco
    Oct 13, 2015 at 17:15
  • This solution may work on my machine, but it is not OS portable.
    – Kalimaha
    Oct 14, 2015 at 6:32
  • @Kalimaha 3 You need a path and back slash as separator at cmd expression. Please, see my edited answer. It works.
    – xunilk
    Oct 14, 2015 at 10:30
  • @xunilk I will try your first solution (and many thanks for it), but not the edit, to avoid the use of os.system which I don't find very elegant.
    – Kalimaha
    Oct 14, 2015 at 11:36
  • @Kalimaha Ok. My first solution really works. I have used it in my own plugins.
    – xunilk
    Oct 14, 2015 at 15:09

1 Answer 1

2

If you are writing a QGIS plugin you can use QgsRasterCalculator and QgsRasterCalculatorEntry classes. For example, this code works well at Python Console (and it could be adapted perfectly for your plugin as a function):

from qgis.core import *
from qgis.gui import *
from qgis.utils import  iface
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

mc = iface.mapCanvas()

Layers = mc.layers()

entries = []

for lyr in Layers:
    # Define bands
    layer = QgsRasterCalculatorEntry()
    layer.ref = lyr.name() +'@1'
    layer.raster = lyr
    layer.bandNumber = 1
    entries.append( layer )

root = '/home/zeito/pyqgis_data/'

expression = '(  ' + entries[1].ref + '   +  '  + entries[0].ref + '  ) / 2 '

print expression

calc = QgsRasterCalculator(expression, 
                            root + 'average.tif', 
                            'GTiff', 
                            Layers[0].extent(), 
                            Layers[0].width(), 
                            Layers[0].height(), 
                            entries )

calc.processCalculation()

I tried it out with two raster (b3_ref and b4_ref); as it can see at the next image. After running the code, the resulting raster was named as average.tif. I corroborated, with the Value Tool Plugin, that its values were as expected. At the Python Console was also printed the calculation expression for purposes of verification.

enter image description here

Editing Note 1:

By using gdal_calc.py, next code works well in the Python Console (path is in Unix notation and the slash, '\', it is necessary at cmd expression).

import os

root_path = '/home/zeito/pyqgis_data/'

cmd = "gdal_calc.py -A " + \
      root_path + \
      "b3_ref.tif -B " + \
      root_path + \
      "b4_ref.tif --outfile=" + \
      root_path + \
      "AVG.tif --cal='(A+B)/2'"

os.system(cmd)

Editing Note 2:

Debugged code. It works in my Python Console. Wrong or unnecessary lines were commented.

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

def average_layers(layers):   #I changed the name
                              #I eliminated "self" to avoid create a plugin

    entries = []

    for l in layers:
#        fileInfo = QFileInfo(l) #external function
#        baseName = fileInfo.baseName() #external function
        layer = QgsRasterCalculatorEntry()
#        QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')
        layer = QgsRasterCalculatorEntry()
        layer.ref = l.name() +'@1'
        layer.raster = l
        layer.bandNumber = 1
        entries.append(layer)

#    for e in entries:
#        QgsMessageLog.logMessage('e: ' + str(e), 'Geobricks TRMM')
#
#    fileInfo = QFileInfo(layers[0]) #external function
#    baseName = fileInfo.baseName() #external function
#    rlayer = QgsRasterLayer(layers[0], baseName) #external function
#    QgsMessageLog.logMessage('rlayer: ' + str(rlayer), 'Geobricks TRMM')

    expression = '(' + entries[1].ref + ' + ' + entries[0].ref + ') / 2'
#    QgsMessageLog.logMessage('Expression: ' + str(expression), 'Geobricks TRMM')
    calc = QgsRasterCalculator(expression,
                               '/home/zeito/pyqgis_data/average.tif', #my absolute pass
                               'GTiff',
                               layers[0].extent(),
                               layers[0].width(),
                               layers[0].height(),
                               entries)

    calc.processCalculation()
#    QgsMessageLog.logMessage('calc: ' + str(calc), 'Geobricks TRMM')

mapcanvas = iface.mapCanvas() 

layers = mapcanvas.layers() #it has the references to the layers
                            #you don't need fileInfo or baseName

average_layers(layers)
0

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.