7

In the QGIS script editor, how do I calculate the main stats (minimum, maximum, mean, median, standard deviation,....) for a single band raster?

I am trying to make a tool for extracting a MASK (0-1) to cover clouds and cloud shadows for Landsat imagery. The algorithm is based on the 'MEAN' and 'Standard Deviation' of Bands BLUE and NIR. To do this I wrote the following script using the 'script editor' within 'Script' tool into the Processing Toolbox of the QGIS. Indeed the script works well, but I would like to automatically calculate 'MEAN' and 'Standard Deviation' instead of entering the values manually.


 ##Shadow-Cloudless MASK=name 
 ##BLUE=raster 
 ##NIR=raster 
 ##MeanBLUE=number 0.0 
 ##StDevBLUE=number 0.0 
 ##MeanNIR=number 0.0 
 ##StDevNIR=number 0.0 
 ##StDevDivisionFactorBLUE=number 1.0 
 ##StDevDivisionFactorNIR=number 1.0 
 ##MASK=output raster

 from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry 
 import numpy as np

 #Get layer object
 layer1 = processing.getObject(BLUE)

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

 #Get number
 number1 = MeanBLUE

 #Get number
 number2 = StDevBLUE

 #Get number
 number3 = MeanNIR

 #Get number
 number4 = StDevNIR

 #Get number
 number5 = StDevDivisionFactorBLUE

 #Get number
 number6 = StDevDivisionFactorNIR

 def mask (BLUE,NIR,MeanBLUE,StDevBLUE,MeanNIR,StDevNIR,StDevDivisionFactorBLUE,StDevDivisionFactorNIR,output):
    entries=[]

    #define raster 1 ("BLUE")
    raster1=QgsRasterCalculatorEntry()
    raster1.ref='BLUE@1'
    raster1.raster=BLUE
    raster1.bandNumber=1
    entries.append(raster1)

    #define raster 2 ("NIR")
    raster2=QgsRasterCalculatorEntry()
    raster2.ref='NIR@1'
    raster2.raster=NIR
    raster2.bandNumber=1
    entries.append(raster2)

    number1='MeanBLUE'
    number2='StDevBLUE'
    number3='StDevDivisionFactorBLUE'
    number4='MeanNIR'
    number5='StDevNIR'
    number6='StDevDivisionFactorNIR'

    #MASK Processing
    calc=QgsRasterCalculator('((("BLUE@1"<('+str(MeanBLUE)+'+('+str(StDevBLUE)+'/'+str(StDevDivisionFactorBLUE)+')))+("NIR@1">('+str(MeanNIR)+'-('+str(StDevNIR)+'/'+str(StDevDivisionFactorNIR)+'))))=' + str(2) + ')',output,'GTiff',NIR.extent(),NIR.width(),NIR.height(),entries)
    calc.processCalculation()

mask(layer1,layer2,number1,number2,number3,number4,number5,number6,MASK)
7

You can use something like the following:

layer = qgis.utils.iface.activeLayer()
provider = layer.dataProvider()
ext = layer.extent()
stats = provider.bandStatistics(1,QgsRasterBandStats.All,ext,0)

print "minimum value = ", stats.minimumValue
print "maximum value = ", stats.maximumValue 
print "mean = ", stats.mean
print "stdDev = ", stats.stdDev

Here it is running from the Script Editor:

Script Editor


EDIT:

I've modified your script and ran it through the Processing Toolbox which worked. It was tested it on a sample raster and an output was generated so hopefully it will work for you:

##Shadow-Cloudless MASK=name 
##BLUE=raster 
##NIR=raster 
##StDevDivisionFactorBLUE=number 1.0 
##StDevDivisionFactorNIR=number 1.0 
##MASK=output raster

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry 
from qgis.core import QgsRasterBandStats
import numpy as np

#Get layer object
layer1 = processing.getObject(BLUE)

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

#Get band statistics of layer1
provider1 = layer1.dataProvider()
ext1 = layer1.extent()
stats1 = provider1.bandStatistics(1,QgsRasterBandStats.All,ext1,0)

#Get band statistics of layer2
provider2 = layer2.dataProvider()
ext2 = layer2.extent()
stats2 = provider2.bandStatistics(1,QgsRasterBandStats.All,ext2,0)

#Get number
number1 = stats1.mean

#Get number
number2 = stats1.stdDev

#Get number
number3 = stats2.mean

#Get number
number4 = stats2.stdDev

#Get number
number5 = StDevDivisionFactorBLUE

#Get number
number6 = StDevDivisionFactorNIR

def mask (BLUE,NIR,MeanBLUE,StDevBLUE,MeanNIR,StDevNIR,StDevDivisionFactorBLUE,StDevDivisionFactorNIR,output):
   entries=[]

   #define raster 1 ("BLUE")
   raster1=QgsRasterCalculatorEntry()
   raster1.ref='BLUE@1'
   raster1.raster=BLUE
   raster1.bandNumber=1
   entries.append(raster1)

   #define raster 2 ("NIR")
   raster2=QgsRasterCalculatorEntry()
   raster2.ref='NIR@1'
   raster2.raster=NIR
   raster2.bandNumber=1
   entries.append(raster2)

   number1='MeanBLUE'
   number2='StDevBLUE'
   number3='StDevDivisionFactorBLUE'
   number4='MeanNIR'
   number5='StDevNIR'
   number6='StDevDivisionFactorNIR'

   #MASK Processing
   calc=QgsRasterCalculator('((("BLUE@1"<('+str(MeanBLUE)+'+('+str(StDevBLUE)+'/'+str(StDevDivisionFactorBLUE)+')))+("NIR@1">('+str(MeanNIR)+'-('+str(StDevNIR)+'/'+str(StDevDivisionFactorNIR)+'))))=' + str(2) + ')',output,'GTiff',NIR.extent(),NIR.width(),NIR.height(),entries)
   calc.processCalculation()

mask(layer1,layer2,number1,number2,number3,number4,number5,number6,MASK)
  • Actually it works perfectly within the python console of the QGIS, but not for the Script Editor, which is what I am looking for. Any solution? – Nico_77 Apr 28 '16 at 13:48
  • @Nico_77 - Huh, works fine for me in the Script Editor. Do you receive any errors? Did you select the raster when you ran it? – Joseph Apr 28 '16 at 13:50
  • it gives me the following message in the log: name 'qgis' is not defined – Nico_77 Apr 28 '16 at 14:17
  • Try adding import qgis at the beginning. – Joseph Apr 28 '16 at 14:19
  • 2
    That is ok. Indeed within the "Style box" into the properties window, if you set "Min/max" and the Accuracy to "Actual(slower)" instead of "Estimate(faster)" the values of the basics stats are the same!! GREAT!! – Nico_77 May 3 '16 at 12:59

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.