3

I'm running a python code (in QGIS 2.18.3) that processes a SPOT image. The code first doing atmospheric correction and then creates NDVI and EVI. The atmospheric correction works great but I can't figure out why I can't create NDVI and EVI. When I run the code, I don't get any errors, the corrected atmospheric bands are created, but that's it, no NDVI or EVI are created. The idea of my code is that after creating the atmospheric corrected bands, the code iterate through the folder and "collects" the necessary atmospheric corrected rasters (as QgsRasterLayer) in order to use them in raster calculator. I've tried to indent the raster calculator syntax in different ways with no success. The if statements at the end works fine, when I use only the first if statement for example, and then use raster calculator it works fine.

I'm attaching here the entire code, but the problem is just in the "#create ndvi, evi and composite bands" section.

##calibrate SPOT 6-7 bands v3=name
##Select_directory=Folder

import glob, os, os.path, math, datetime
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsRasterLayer, QgsRectangle, QgsVectorLayer

date_format = '%Y-%m-%d'
TIME = "<TIME>"
GAIN="GAIN"
BIAS="BIAS"
solar_irradiance="Solar irradiance"
SUN_ELEVATION = "SUN_ELEVATION"
band_name = ["Red", "Green", "Blue", "NIR"]

for file in os.listdir(Select_directory):
    if file[:3]=="DIM":

        xml_file = os.path.join(Select_directory,file) 
        z = open(xml_file,'r')
        Sun_elevation_list=[]
        Gain_list =[]
        Bias_list=[]
        solar_irradiance_value_list = []
        solar_irradiance_line_list = []

        for line, text in enumerate(z):  
            if GAIN in text:   
                y = text.split(">")
                yy = y[1].split("<")
                if len((yy[0]))<6:
                    Gain = yy[0]
                    Gain_list.append(Gain)

            if BIAS in text:   
                y = text.split(">")
                yy = y[1].split("<")
                if len((yy[0]))<6:
                    Bias= yy[0]
                    Bias_list.append(Bias)

            if solar_irradiance in text:
                solar_irradiance_line_list.append(line+3) # the plus 3 is because the actual value is 3 rows below 

            if SUN_ELEVATION in text:
                s = text.split("<SUN_ELEVATION>")
                ss = s[1].split("<")
                sss = float(ss[0])
                Sun_elevation_list.append(sss)

            if TIME in text:
                t = text.split("<TIME>")
                date = t[1][:10]
                dt = datetime.datetime.strptime(date, date_format)
                tt = dt.timetuple()
                Julian_Day = tt.tm_yday

        f=open(xml_file)
        lines=f.readlines()
        for t in range(len(solar_irradiance_line_list)):
            text1 = lines[solar_irradiance_line_list[t]]
            v =text1.split("VALUE>")
            vv = v[1].split("<")
            solar_irradiance_value_list.append(vv[0])

        distance = 1-0.01672*math.cos(math.radians(0.9856*(Julian_Day-4)))
        distance_square = distance*distance
        SUN_ELEVATION_AVG = sum(Sun_elevation_list)/len(Sun_elevation_list)
        solar_zenithal_angle_celsius = 90-SUN_ELEVATION_AVG
        solar_zenithal_angle_radian = math.radians(solar_zenithal_angle_celsius)
        solar_cos = math.cos(solar_zenithal_angle_radian)
        pi = math.pi


#create TOA seperate bands
os.chdir(Select_directory)
for raster in glob.glob("*.tif"):
        if raster[:3]=="IMG":
            fileInfo = QFileInfo(raster)
            baseName = fileInfo.baseName()
            rlayer = QgsRasterLayer(raster, baseName)
            band_number = 1

            for n in range(len(Gain_list)):
                syntax = '10000*({}*{}*((A/{})+{}))/({}*{})'.format(pi,distance_square,Gain_list[n],Bias_list[n],solar_irradiance_value_list[n],solar_cos) 
                output_raster = os.path.join(Select_directory,"band_"+band_name[n])
                outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', rlayer,band_number,None,'1',None,'1',None,'1',None,'1',None,'1',syntax,'0',4,None,output_raster)
                band_number+=1

#create ndvi, evi and composite bands
os.chdir(Select_directory)
NDVI = os.path.join(Select_directory, "NDVI")
EVI = os.path.join(Select_directory, "EVI")
NDVI_syntax='(A-B)/(A+B)'
EVI_syntax= '2.5*(A-B)/(A+6*B-7.5*C+1)'

for raster in glob.glob("*.tif"):
    fileInfo = QFileInfo(raster)
    baseName = fileInfo.baseName()

    if baseName == "band_Red":
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        RED = QgsRasterLayer(raster, baseName)
    elif baseName == "band_Green":
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        GREEN = QgsRasterLayer(raster, baseName)
    elif baseName == "band_Blue":
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        BLUE = QgsRasterLayer(raster, baseName)
    elif baseName == "band_NIR":
        fileInfo = QFileInfo(raster)
        baseName = fileInfo.baseName()
        NIR = QgsRasterLayer(raster, baseName)


outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', NIR,'1',RED,'1',None,'1',None,'1',None,'1',None,'1',None,'1',NDVI_syntax,'0',4,None,NDVI) 
outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', NIR,'1',RED,'1',BLUE,'1',None,'1',None,'1',None,'1',None,'1',None,'1',EVI_syntax,'0',4,None,EVI)  
2
  • If your problem is "just in [some] section" then reduce your code to create a minimal example. stackoverflow.com/help/mcve This will potentially lead you to the solution already. Mar 6, 2017 at 9:52
  • You have more parameters to be used in your command lines (some of them in a wrong position). Please, see my answer.
    – xunilk
    Mar 6, 2017 at 10:23

1 Answer 1

2

This is the complete list of parameters to be used in processing 'gdalogr:rastercalculator'.

ALGORITHM: Raster calculator
    INPUT_A <ParameterRaster>
    BAND_A <ParameterString>
    INPUT_B <ParameterRaster>
    BAND_B <ParameterString>
    INPUT_C <ParameterRaster>
    BAND_C <ParameterString>
    INPUT_D <ParameterRaster>
    BAND_D <ParameterString>
    INPUT_E <ParameterRaster>
    BAND_E <ParameterString>
    INPUT_F <ParameterRaster>
    BAND_F <ParameterString>
    FORMULA <ParameterString>
    NO_DATA <ParameterString>
    RTYPE <ParameterSelection>
    EXTRA <ParameterString>
    OUTPUT <OutputRaster>


RTYPE(Output raster type)
    0 - Byte
    1 - Int16
    2 - UInt16
    3 - UInt32
    4 - Int32
    5 - Float32
    6 - Float64

So, your code lines:

outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', NIR,'1',RED,'1',None,'1',None,'1',None,'1',None,'1',None,'1',NDVI_syntax,'0',4,None,NDVI) 
outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', NIR,'1',RED,'1',BLUE,'1',None,'1',None,'1',None,'1',None,'1',None,'1',EVI_syntax,'0',4,None,EVI)

are wrong (and you didn't import processing python module).

To test my approach I used next code employing your terminology and with atmospheric corrected NIR and RED bands:

import processing

NIR = processing.getObjectFromName('NIR')
RED = processing.getObjectFromName('RED')

NDVI_syntax = '(A-B)/(A+B)'

outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', 
                                                    NIR,         #INPUT_A <ParameterRaster>
                                                    '1',         #BAND_A <ParameterString>
                                                    RED,         #INPUT_B <ParameterRaster>
                                                    '1',         #BAND_B <ParameterString>
                                                    None,        #INPUT_C <ParameterRaster>
                                                    '1',         #BAND_C <ParameterString>
                                                    None,        #INPUT_D <ParameterRaster>
                                                    '1',         #BAND_D <ParameterString>
                                                    None,        #INPUT_E <ParameterRaster>
                                                    '1',         #BAND_E <ParameterString>
                                                    None,        #INPUT_F <ParameterRaster>
                                                    '1',         #BAND_F <ParameterString>
                                                    NDVI_syntax, #FORMULA <ParameterString>
                                                    '',          #NO_DATA <ParameterString>
                                                    5,           #RTYPE <ParameterSelection>
                                                    '0',         #EXTRA <ParameterString>
                                                    None)        #OUTPUT <OutputRaster>

NDVI = QgsRasterLayer(outputs_GDALOGRRASTERCALCULATOR_1['OUTPUT'],
                      'NDVI')

QgsMapLayerRegistry.instance().addMapLayer(NDVI)

Afterward running the code at the Python Console of QGIS, NDVI layer was loaded at the Map Canvas:

enter image description here

NDVI values (between -1 and 1) were classified with the ranges suggested in this link: https://en.wikipedia.org/wiki/Normalized_Difference_Vegetation_Index

Blue area corresponds effectively to a lake. It worked correctly.

2
  • Thank you very much, I've deleted the redundant none-s in my code and it works. When running python code from QGIS script editor it works without the need to specifically import processing
    – user88484
    Mar 6, 2017 at 15:20
  • @xunilk How can I get those NDVI values (-1,0.3,+1 in above image) using python code? Apr 9, 2017 at 6:11

Your Answer

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

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