QGIS secuencia de comandos para calcular estadísticas para una sola banda de ráster

Question

Alguien me puede ayudar a escribir un código de secuencia de comandos, en QGIS editor de secuencias de comandos, para el cálculo de las principales estadísticas (mínimo, máximo, media, mediana, desviación estándar,....) para una sola banda de trama?

En realidad estoy tratando de hacer una herramienta para la extracción de una MÁSCARA (0-1) a la cubierta de nubes y la nube de sombras para las imágenes Landsat. El algoritmo se basa en el "MEDIO" y "Desviación Estándar" de las Bandas AZUL y NIR. Para ello he escrito la siguiente secuencia de comandos utilizando el "editor de secuencias de comandos' dentro de 'Script' de la herramienta en el Procesamiento de la caja de herramientas de QGIS. De hecho, el script funciona bien, pero me gustaría para calcular automáticamente el "medio" y "Desviación Estándar" en lugar de introducir los valores manualmente. Alguna idea?

---

 ##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)

Answer 1

Usted puede usar algo como la siguiente:

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

---

Aquí se está ejecutando desde el Editor de secuencias de Comandos:

---

EDITAR:

He modificado el script y corrió a través del Procesamiento de la caja de herramientas que trabajó. Fue probado en una muestra de trama y una salida se generó así que espero que va a trabajar para usted:

##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)