[GIS] Updating contrast enhancements for raster layers using bandStatistics in PyQGIS

pyqgisrasterstretch

Found this function to update statistics for a raster layer :

rlayer.setContrastEnhancement(QgsContrastEnhancement.StretchToMinimumMaximum,
                          QgsRaster.ContrastEnhancementMinMax)

Which works great for minmax and cumulative but stddev is hardcoded to 1 standard deviation. For cumulative it reads the lower and upper values from the settings but no setting for stddev.

So I have looked at the C++ code in the api and tried to convert that into a QGIS script.
The only problem is that the band statistics and cumulative cut functions are returning "nan". Take a look at my code segment below and see if you can help:

       myEnhancement = QgsContrastEnhancement(myType)
       myEnhancement.setContrastEnhancementAlgorithm(ContrastEnhancement,True)

       if Limits == 0:
          myRasterBandStats = rlayer.renderer().bandStatistics( myBand,
                                     QgsRasterBandStats.Min | QgsRasterBandStats.Max)
          myMin = myRasterBandStats.minimumValue
          myMax = myRasterBandStats.maximumValue

       elif Limits == 1:
            myRasterBandStats = rlayer.renderer().bandStatistics( myBand,
                                   QgsRasterBandStats.Mean | QgsRasterBandStats.StdDev)
            myMin = myRasterBandStats.mean - ( StdDev * myRasterBandStats.stdDev )
            myMax = myRasterBandStats.mean + ( StdDev * myRasterBandStats.stdDev )

       elif Limits == 2:  
            myMin,myMax = rlayer.renderer().cumulativeCut( myBand, CumulativeLower,
                                                          CumulativeUpper)

       raise GeoAlgorithmExecutionException( str(myBand) + " " + str(myMin) + " " + str(myMax))
       myEnhancement.setMinimumValue( myMin )
       myEnhancement.setMaximumValue( myMax )

       if renderType == "singlebandgray":
          rlayer.renderer().setContrastEnhancement(myEnhancement)
       elif renderType == "multibandcolor":
          if Band == 0:
             rlayer.renderer().setRedContrastEnhancement(myEnhancement)
          elif Band == 1:
             rlayer.renderer().setGreenContrastEnhancement(myEnhancement)
          elif Band == 2:
             rlayer.renderer().setBlueContrastEnhancement(myEnhancement)
       rlayer.triggerRepaint()

Heres the full script using the first function which works great but does not enable you to have more than 1 standard deviation

##[Example scripts]=group
##Limits=selection "MinMax";"StdDev";"Cumulative"
##Stretch=selection "NoStretch";"StretchToMinMax";"StretchAndClipToMinMax";"ClipToMinMax"

import os
from qgis.core import *
from qgis.gui import *
import qgis.utils

layers = qgis.utils.iface.legendInterface().layers()

for layer in layers:
    layerType = layer.type()
    if layerType == QgsMapLayer.RasterLayer:
       if Stretch==0:
          ContrastEnhancement = QgsContrastEnhancement.NoEnhancement
       elif Stretch==1:
            ContrastEnhancement = QgsContrastEnhancement.StretchToMinimumMaximum
       elif Stretch==2:
            ContrastEnhancement = QgsContrastEnhancement.StretchAndClipToMinimumMaximum
       elif Stretch==3:
            ContrastEnhancement = QgsContrastEnhancement.ClipToMinimumMaximum

       if Limits==0:
          layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementMinMax)
       elif Limits==1:
          layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementStdDev)
       elif Limits==2:
           layer.setContrastEnhancement(ContrastEnhancement,QgsRaster.ContrastEnhancementCumulativeCut)
       layer.triggerRepaint()

Best Answer

I fixed the script, you need to call bandStatistics on the dataProvider() not the renderer().So heres the final solution that works well

##[Example scripts]=group
##Limits=selection "MinMax";"StdDev";"Cumulative"
##Stretch=selection "NoStretch";"StretchToMinMax";"StretchAndClipToMinMax";"ClipToMinMax"
##StdDev=number 1.0
##CumulativeLower=number 0.02
##CumulativeUpper=number 0.98

import os
from qgis.core import *
from qgis.gui import *
import qgis.utils

from processing.core.GeoAlgorithmExecutionException import \
        GeoAlgorithmExecutionException

layers = qgis.utils.iface.legendInterface().layers()

for layer in layers:
    layerType = layer.type()
    if layerType == QgsMapLayer.RasterLayer:
       renderType = layer.renderer().type()
       if renderType == "multibandcolor":
          bands = 3
       else : bands = 1

       for Band in range(bands):
           if renderType == "multibandcolor":
              if Band == 0:
                 myBand = layer.renderer().redBand()
              elif Band == 1:
                 myBand = layer.renderer().greenBand()
              elif Band == 2:
                 myBand = layer.renderer().blueBand()
           else : myBand = layer.renderer().grayBand()
           myType = layer.renderer().dataType(myBand)

           if Stretch==0:
              ContrastEnhancement = QgsContrastEnhancement.NoEnhancement
           elif Stretch==1:
              ContrastEnhancement = QgsContrastEnhancement.StretchToMinimumMaximum
           elif Stretch==2:
             ContrastEnhancement = QgsContrastEnhancement.StretchAndClipToMinimumMaximum
           elif Stretch==3:
             ContrastEnhancement = QgsContrastEnhancement.ClipToMinimumMaximum

           myEnhancement = QgsContrastEnhancement(myType)
           myEnhancement.setContrastEnhancementAlgorithm(ContrastEnhancement,True)

           if Limits == 0:
              myRasterBandStats = layer.dataProvider().bandStatistics( myBand, QgsRasterBandStats.Min | QgsRasterBandStats.Max)
              myMin = myRasterBandStats.minimumValue
              myMax = myRasterBandStats.maximumValue

           elif Limits == 1:
                myRasterBandStats = layer.dataProvider().bandStatistics( myBand, QgsRasterBandStats.Mean | QgsRasterBandStats.StdDev)
                myMin = myRasterBandStats.mean - ( StdDev * myRasterBandStats.stdDev )
                myMax = myRasterBandStats.mean + ( StdDev * myRasterBandStats.stdDev )

           elif Limits == 2:  
                myMin,myMax = layer.dataProvider().cumulativeCut( myBand, CumulativeLower, CumulativeUpper)

#           raise GeoAlgorithmExecutionException( str(myBand) + " " + str(myMin) + " " + str(myMax))
           myEnhancement.setMinimumValue( myMin )
           myEnhancement.setMaximumValue( myMax )


           if renderType == "multibandcolor":
              if Band == 0:
                 layer.renderer().setRedContrastEnhancement(myEnhancement)
              elif Band == 1:
                 layer.renderer().setGreenContrastEnhancement(myEnhancement)
              elif Band == 2:
                 layer.renderer().setBlueContrastEnhancement(myEnhancement)
           else : layer.renderer().setContrastEnhancement(myEnhancement)
       layer.triggerRepaint()
Related Question