3

I'm trying to write a Python script in QGIS using the Python console that applies the NDVI formula to all rasters in the Layers Panel, and outputs a new NDVI raster for each original raster. My skill level in Python is very minimal, considering that, I've made some decent progress thus far. My current script references all rasters in the Layers Panel, and outputs a resulting raster with a pre-determined naming convention.

The problem that I'm having is the new NDVI rasters have "nan" values and are therefore useless. I've read a number of blogs on this issue but to no avail. I believe the issue exists in how I'm referencing the red and nir bands (see the Raster Calculator and how the resulting NDVI layer only has one band), and/or how they're being applied in the Raster Calculator.

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

    lddLrs = qgis.utils.iface.legendInterface().layers()
    path = "D:/QGIS/Project/"
    for lyr in lddLrs:
        entries = []
        ras = QgsRasterCalculatorEntry()
        ras.ref = 'ras@1'
        ras.raster = lyr
        ras.bandNumber = 1
        ras = QgsRasterCalculatorEntry()
        ras.ref = 'ras@4'
        ras.raster = lyr
        ras.bandNumber = 4
        entries.append( ras )
        calc = QgsRasterCalculator( '(ras@4 -  ras@1) / (ras@4 + ras1)', path 
        + lyr.name() + "_NDVI.tif", 'GTiff', lyr.extent(), lyr.width(), 
        lyr.height(), entries )
        calc.processCalculation()

enter image description here

1 Answer 1

1

So after some experimenting, I finally got the code to work as needed.

    lddLrs = qgis.utils.iface.legendInterface().layers()
    path = "D:/QGIS/Project/"
    for lyr in lddLrs:
        entries = []
        #Define red band#
        ras1 = QgsRasterCalculatorEntry()
        ras1.ref = 'layer.name() +@1'
        ras1.raster = lyr
        ras1.bandNumber = 1
        entries.append(ras1)
        #Define nir band#
        ras2 = QgsRasterCalculatorEntry()
        ras2.ref = 'layer.name() +@4'
        ras2.raster = lyr
        ras2.bandNumber = 4
        entries.append( ras2 )
        calc = QgsRasterCalculator( '("layer.name() +@4" -  "layer.name() 
       +@1") / ("layer.name() +@4" + "layer.name() +@1")', path + lyr.name() 
       + "_NDVI.tif", 'GTiff', lyr.extent(), lyr.width(), lyr.height(), 
       entries )
       calc.processCalculation()

Boom.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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