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)
