# Uses Landsat data as downloaded from earthexplorer.usgs.gov
# The script automatically handles Landsat 5, 7 and 8
# Files are extracted, radiation and toa relfectance is calculated as described in the repsective Landsat Handbooks

import tarfile
import os, shutil
import math
import numpy as np


# The following options allow control over the desired output and handling of intermediate data

deleteOriginal = 3 # 1 = keep all, 2 = delete all, 3 = delete but keep metafile
# Note: folder needs to be empty at start, else all is deleted
saveRadiation = False # if True Radiation Tiffs will be saved
saveReflectance = 1 # if True Reflectance Tiffs will be saved
calcKelvin = 0 # 0 = no calculation, 1 = calc Band 10&11 or 61/62,
# 2 = calc B10/61 only, 3=calc B11/62 only
allIndices = 1 # 0 = none, 1= all, 2 = NDVI, 3 = EVI, 4 = SAVI,
# 5 = MSAVI, 6 = NDMI, 7 = NBR, 8 = NBR2


############################################################################
#All of the follwing are used to calculate earth-sun distance, adapted from
#ttp://davit.ece.vt.edu/davitpy/_modules/utils/calcSun.html#calcSunRadVector
def calcGeomMeanAnomalySun( t ):
""" Calculate the Geometric Mean Anomaly of the Sun (in degrees) """
M = 357.52911 + t * ( 35999.05029 - 0.0001537 * t)
return M # in degrees

def calcEccentricityEarthOrbit( t ):
""" Calculate the eccentricity of earth's orbit (unitless) """
e = 0.016708634 - t * ( 0.000042037 + 0.0000001267 * t)
return e # unitless

def calcSunEqOfCenter( t ):
"""Calculate the equation of center for the sun (in degrees) """
mrad = np.radians(calcGeomMeanAnomalySun(t))
sinm = np.sin(mrad)
sin2m = np.sin(mrad+mrad)
sin3m = np.sin(mrad+mrad+mrad)
C = (sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) +
sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289)
return C # in degrees

def calcSunTrueAnomaly( t ):
m = calcGeomMeanAnomalySun(t)
c = calcSunEqOfCenter(t)
v = m + c
return v # in degrees

# Calculate earth-sun distance in AE, adapted from
def ES_dist(t):

#eccent = 0.01672592 # Eccentricity of Earth orbit
#axsmaj = 1.4957 # Semi-major axis of Earth orbit (km)
#solyr = 365.2563855 # Number of days in a solar year

v = calcSunTrueAnomaly(t)
e = calcEccentricityEarthOrbit(t)
R = (1.000001018 * (1. - e * e)) / ( 1. + e * np.cos( np.radians(v) ) )
return R
############################################################################


# Define functions for spectral indices
def NDVI(R,NIR):
#LS8 = B4,B5 LS57 = B3,B4
ndvi = (NIR - R) / (NIR + R)
ndvi = ndvi.astype(np.float32)
return ndvi

def EVI(B,R,NIR,L=1):
#LS8 = B2,B4,B5 LS57 = B1,B3,B4
evi = (NIR - R) / (NIR + 6*R - 7.5 * B + L)
evi = evi.astype(np.float32)
return evi

def SAVI(R,NIR,L=0.5):
#LS8 = B4,B5 LS57 = B3,B4
savi = ((NIR-R) / (NIR + R + L)) * (1+L)
savi = savi.astype(np.float32)
return savi

def MSAVI(R,NIR):
#LS8 = B4,B5 LS57 = B3,B4
msavi = (2 * NIR + 1 - np.sqrt((2*NIR+1)**2 - 8 * (NIR - R))) / 2
msavi = msavi.astype(np.float32)
return msavi

def NDMI(NIR,SWIR):
#LS8 = B5,B6 LS57 = B4,B5
ndmi = (NIR-SWIR) / (NIR+SWIR)
ndmi = ndmi.astype(np.float32)
return ndmi

def NBR(NIR,SWIR):
#LS8 = B5,B7 LS57 = B4,B7
nbr = (NIR - SWIR) / (NIR+SWIR)
nbr = nbr.astype(np.float32)
return nbr

def NBR2(SWIR1,SWIR2):
#LS8 = B6,B7 LS57 = B5,B7
nbr2 = (SWIR1 - SWIR2) / (SWIR1 + SWIR2)
nbr2 = nbr2.astype(np.float32)
return nbr2


# function checks if curDir already exists -> does nothing, else creates it
def chkdir2(curDir):
if os.path.isdir(curDir):
return
else:
try:
os.mkdir(curDir)
except:
print("Cannot create Folder: ", curDir)

return curDir



# function reads one raster (path) as an array and returns array object
# input format is GeoTiff
def singleTifToArray(inRas):

firstRasGDAL = gdal.Open(inRas, GA_ReadOnly)
cols = firstRasGDAL.RasterXSize
rows = firstRasGDAL.RasterYSize

dataset = gdal.Open(inRas, GA_ReadOnly)
array = dataset.ReadAsArray(0, 0, cols, rows)

return array



# Function that converts a numpy array to a GeoTIFF and saves it to disk, needs an existing GeoTIFF file as a blueprint
# Changed after the original
# http://gis.stackexchange.com/questions/58517/python-gdal-save-array-as-raster-with-projection-from-other-file
# inTiff is an exisiting GeoTiff file, the attributes from this file are used to create the output
# array is the array that will be saved as a tiff
# outFile is the path and name of the desired output GeoTiff

def array_to_raster(inTiff,array,outFile):

inDataset = gdal.Open(inTiff, GA_ReadOnly)

# You need to get those values like you did.
x_pixels = inDataset.RasterXSize # number of pixels in x
y_pixels = inDataset.RasterYSize # number of pixels in y
PIXEL_SIZE = inDataset.GetGeoTransform()[1] # size of the pixel...
x_min = inDataset.GetGeoTransform()[0]
y_max = inDataset.GetGeoTransform()[3] # x_min & y_max are like the "top left" corner.
wkt_projection = inDataset.GetProjectionRef()

driver = gdal.GetDriverByName('GTiff')

outDataset = driver.Create(
outFile,
x_pixels,
y_pixels,
1,
gdal.GDT_Float32, )

outDataset.SetGeoTransform((
x_min, # 0
PIXEL_SIZE, # 1
0, # 2
y_max, # 3
0, # 4
-PIXEL_SIZE))

outDataset.SetProjection(wkt_projection)
outDataset.GetRasterBand(1).WriteArray(array)
outDataset.FlushCache() # Write to disk.
return outDataset, outDataset.GetRasterBand(1)


# Acutal funtion that will calculate the output rasters, inFile must be an original *.tar.gz file
def Landsat(inFile,outFol,deleteOriginal,saveRadiation,saveReflectance,calcKelvin,allIndices):

#Extract Metainfo from packaged Landsat Data
tar = tarfile.open(inFile)
members = tar.getmembers()

#Name of current Landsat file, extracted from first band name
LSname = inFile[inFile.rfind("/")+1:inFile.find(".")]
julDay = int(LSname[13:16])
year = int(LSname[9:13])

#Define output Folder
curFol = outFol + LSname + "/"
chkdir2(curFol)

tar.extractall(curFol)

#collect filenames in a list to maybe delete later
orgList = []
for x in os.listdir(curFol):
orgList.append(curFol + x)


#Handle Landsat 8 differently than 5 or 7
if LSname[2] == '8':

#Find the metafile and extract variables
for x in os.listdir(curFol):
#Following will only consider Landsat8 as valid input
#http://landsat.usgs.gov/Landsat8_Using_Product.php
if x[-7:].lower() == "mtl.txt":
#metaFile = tar.extractfile(x) #metaFiel is now a fileObject
metaFile = open(curFol + x, "r")

for lines in metaFile.readlines():
#lines still contain whitespaces
#strLine = lines.decode("utf-8").replace(" ", "")
strLine = lines.replace(" ", "")

if strLine[:13] == 'SPACECRAFT_ID':
spacecraft_ID = strLine[15:-2]
if strLine[:8] == 'WRS_PATH':
wrs_path = strLine[9:]
if strLine[:7] == 'WRS_ROW':
wrs_row = strLine[8:]
if strLine[:11] == 'CLOUD_COVER':
cloud_cover = float(strLine[12:-1])
if strLine[:13] == 'DATE_ACQUIRED':
date_aquired_str = strLine[14:-1]
if strLine[:17] == 'SCENE_CENTER_TIME':
scene_center_time_str = strLine[18:-1]
if strLine[:8] == 'UTM_ZONE':
utm_zone = strLine[9:-1]

if strLine[:11] == 'SUN_AZIMUTH':
sun_azimuth = float(strLine[12:])
if strLine[:13] == 'SUN_ELEVATION':
sun_elevation = float(strLine[14:])
if strLine[:18] == 'EARTH_SUN_DISTANCE':
earth_sun_distance = float(strLine[19:])

#Radiances, first ones need to make sure Band_10 or Band_11 are not chosen instead
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_1' and strLine[23] == "=":
radiance_max_B1 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_1' and strLine[23] == "=":
radiance_min_B1 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_2':
radiance_max_B2 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_2':
radiance_min_B2 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_3':
radiance_max_B3 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_3':
radiance_min_B3 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_4':
radiance_max_B4 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_4':
radiance_min_B4 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_5':
radiance_max_B5 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_5':
radiance_min_B5 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_6':
radiance_max_B6 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_6':
radiance_min_B6 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_7':
radiance_max_B7 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_7':
radiance_min_B7 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_8':
radiance_max_B8 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_8':
radiance_min_B8 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_9':
radiance_max_B9 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_9':
radiance_min_B9 = float(strLine[24:-2])

if strLine[:24] == 'RADIANCE_MAXIMUM_BAND_10':
radiance_max_B10 = float(strLine[25:-2])
if strLine[:24] == 'RADIANCE_MINIMUM_BAND_10':
radiance_min_B10 = float(strLine[25:-2])
if strLine[:24] == 'RADIANCE_MAXIMUM_BAND_11':
radiance_max_B11 = float(strLine[25:-2])
if strLine[:24] == 'RADIANCE_MINIMUM_BAND_11':
radiance_min_B11 = float(strLine[25:-2])


#Reflectances, only 9 bands in metafile
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_1':
reflectance_max_B1 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_1':
reflectance_min_B1 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_2':
reflectance_max_B2 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_2':
reflectance_min_B2 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_3':
reflectance_max_B3 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_3':
reflectance_min_B3 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_4':
reflectance_max_B4 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_4':
reflectance_min_B4 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_5':
reflectance_max_B5 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_5':
reflectance_min_B5 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_6':
reflectance_max_B6 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_6':
reflectance_min_B6 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_7':
reflectance_max_B7 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_7':
reflectance_min_B7 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_8':
reflectance_max_B8 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_8':
reflectance_min_B8 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MAXIMUM_BAND_9':
reflectance_max_B9 = float(strLine[27:-2])
if strLine[:26] == 'REFLECTANCE_MINIMUM_BAND_9':
reflectance_min_B9 = float(strLine[27:-2])

# Min_Max pixel value
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_1' and strLine[23] == "=":
quant_max_B1 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_1' and strLine[23] == "=":
quant_min_B1 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_2':
quant_max_B2 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_2':
quant_min_B2 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_3':
quant_max_B3 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_3' :
quant_min_B3 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_4' :
quant_max_B4 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_4' :
quant_min_B4 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_5':
quant_max_B5 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_5' :
quant_min_B5 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_6':
quant_max_B6 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_6' :
quant_min_B6 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_7':
quant_max_B7 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_7':
quant_min_B7 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_8' :
quant_max_B8 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_8' :
quant_min_B8 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_9' :
quant_max_B9 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_9':
quant_min_B9 = int(strLine[24:])

if strLine[:24] == 'QUANTIZE_CAL_MAX_BAND_10':
quant_max_B10 = int(strLine[25:])
if strLine[:24] == 'QUANTIZE_CAL_MIN_BAND_10':
quant_min_B10 = int(strLine[25:])
if strLine[:24] == 'QUANTIZE_CAL_MAX_BAND_11':
quant_max_B11 = int(strLine[25:])
if strLine[:24] == 'QUANTIZE_CAL_MIN_BAND_11':
quant_min_B11 = int(strLine[25:])


# Radiometric Rescaling
if strLine[:20] == 'RADIANCE_MULT_BAND_1' and strLine[20] == "=":
radiance_mult_B1 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_2':
radiance_mult_B2 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_3':
radiance_mult_B3 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_4':
radiance_mult_B4 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_5':
radiance_mult_B5 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_6':
radiance_mult_B6 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_7':
radiance_mult_B7 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_8':
radiance_mult_B8 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_9':
radiance_mult_B9 = float(strLine[21:])
if strLine[:21] == 'RADIANCE_MULT_BAND_10':
radiance_mult_B10 = float(strLine[22:])
if strLine[:21] == 'RADIANCE_MULT_BAND_11':
radiance_mult_B11 = float(strLine[22:])

if strLine[:19] == 'RADIANCE_ADD_BAND_1' and strLine[19] == "=":
radiance_add_B1 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_2':
radiance_add_B2 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_3':
radiance_add_B3 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_4':
radiance_add_B4 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_5':
radiance_add_B5 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_6':
radiance_add_B6 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_7':
radiance_add_B7 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_8':
radiance_add_B8 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_9':
radiance_add_B9 = float(strLine[20:])
if strLine[:20] == 'RADIANCE_ADD_BAND_10':
radiance_add_B10 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_ADD_BAND_11':
radiance_add_B11 = float(strLine[21:])

if strLine[:23] == 'REFLECTANCE_MULT_BAND_1' and strLine[23] == "=":
reflectance_mult_B1 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_2':
reflectance_mult_B2 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_3':
reflectance_mult_B3 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_4':
reflectance_mult_B4 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_5':
reflectance_mult_B5 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_6':
reflectance_mult_B6 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_7':
reflectance_mult_B7 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_8':
reflectance_mult_B8 = float(strLine[24:])
if strLine[:23] == 'REFLECTANCE_MULT_BAND_9':
reflectance_mult_B9 = float(strLine[24:])

if strLine[:22] == 'REFLECTANCE_ADD_BAND_1' and strLine[22] == "=":
reflectance_add_B1 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_2':
reflectance_add_B2 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_3':
reflectance_add_B3 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_4':
reflectance_add_B4 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_5':
reflectance_add_B5 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_6':
reflectance_add_B6 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_7':
reflectance_add_B7 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_8':
reflectance_add_B8 = float(strLine[23:])
if strLine[:22] == 'REFLECTANCE_ADD_BAND_9':
reflectance_add_B9 = float(strLine[23:])

#Thermal Constants
if strLine[:19] == 'K1_CONSTANT_BAND_10':
k1_const_B10 = float(strLine[20:])
if strLine[:19] == 'K1_CONSTANT_BAND_11':
k1_const_B11 = float(strLine[20:])
if strLine[:19] == 'K2_CONSTANT_BAND_10':
k2_const_B10 = float(strLine[20:])
if strLine[:19] == 'K2_CONSTANT_BAND_11':
k2_const_B11 = float(strLine[20:])

metaFile.close()
members = []


#calculate radiation and toa relfectance
for x in range(1,10):
#read bands as arrays
exec( 'arrayDN{0} = singleTifToArray(curFol + LSname + "_B{0}.TIF")'.format(x) )

#outStr = curFol + LSname + "_B{0}.TIF".format(x)
#print(outStr)

#convert to radiance and convert to 32-bit floating point for memory saving
exec( 'lambda{0} = radiance_mult_B{0} * arrayDN{0} + radiance_add_B{0}'.format(x) )
exec( 'lambda{0} = lambda{0}.astype(np.float32)'.format(x))

#convert to reflectance and convert to 32-bit floating point for memory saving
exec( 'reflectance{0} = ( (reflectance_mult_B{0} * arrayDN{0} + reflectance_add_B{0}) / \
math.sin(math.radians(sun_elevation)) )'.format(x) )
exec( 'reflectance{0} = reflectance{0}.astype(np.float32)'.format(x))

#empty first array for memory saving
exec( 'arrayDN{0} = []'.format(x) )


#calculate temperature in Kelvin from bands 10 and 11
if calcKelvin > 0:
if calcKelvin == 1 or calcKelvin == 2:
arrayDN10 = singleTifToArray(curFol + LSname + "_B10.TIF")
lambda10 = radiance_mult_B10 * arrayDN10 + radiance_add_B10
lambda10 = lambda10.astype(np.float32)
t10 = k2_const_B10 / ( np.log( (k1_const_B10 / lambda10) + 1) ) #T in Kelvin
t10 = t10.astype(np.float32)

array_to_raster(curFol + LSname + "_B10.TIF", t10,
curFol + "Temp_B10.TIF")
arrayDN10 = []

if calcKelvin == 1 or calcKelvin == 3:
arrayDN11 = singleTifToArray(curFol + LSname + "_B11.TIF")
lambda11 = radiance_mult_B11 * arrayDN11 + radiance_add_B11
lambda11 = lambda11.astype(np.float32)
t11 = k2_const_B11 / ( np.log( (k1_const_B11 / lambda11) + 1) )
t11 = t11.astype(np.float32)

array_to_raster(curFol + LSname + "_B11.TIF", t11,
curFol + "Temp_B11.TIF")
arrayDN11 = []



#Save radiation rasters to disk, arrays will be deleted either way
if saveRadiation:
if calcKelvin == False:
endRange = 10
else:
endRange = 12
for x in range(1,endRange):
exec( 'array_to_raster(curFol + LSname + "_B{0}.TIF", lambda{0}, \
curFol + "Radiation_B{0}.TIF")'.format(x) )
exec( 'lambda{0} = []'.format(x) )
else:
if calcKelvin == False:
endRange = 10
else:
endRange = 12
for x in range(1,endRange):
exec( 'lambda{0} = []'.format(x) )


#Calculate spectral indices
if allIndices > 0:
if allIndices == 1 or allIndices == 2:
ndviAr = NDVI(reflectance4,reflectance5)
array_to_raster(curFol + LSname + "_B4.TIF", ndviAr,
curFol + "Ind_NDVI_" + LSname[:-5] + ".TIF")
ndviAr = []

if allIndices == 1 or allIndices == 3:
eviAr = EVI(reflectance2,reflectance4,reflectance5,L=1)
array_to_raster(curFol + LSname + "_B2.TIF", eviAr,
curFol + "Ind_EVI_" + LSname[:-5] + ".TIF")
eviAr = []

if allIndices == 1 or allIndices == 4:
saviAr = SAVI(reflectance4,reflectance5,L=0.5)
array_to_raster(curFol + LSname + "_B4.TIF", saviAr,
curFol + "Ind_SAVI_" + LSname[:-5] + ".TIF")
saviAr = []

if allIndices == 1 or allIndices == 5:
msaviAr = MSAVI(reflectance4,reflectance5)
array_to_raster(curFol + LSname + "_B4.TIF", msaviAr,
curFol + "Ind_MSAVI_" + LSname[:-5] + ".TIF")
msaviAr = []

if allIndices == 1 or allIndices == 6:
ndmiAr = NDMI(reflectance5,reflectance6)
array_to_raster(curFol + LSname + "_B5.TIF", ndmiAr,
curFol + "Ind_NDMI_" + LSname[:-5] + ".TIF")
ndmiAr = []

if allIndices == 1 or allIndices == 7:
nbrAr = NBR(reflectance5,reflectance7)
array_to_raster(curFol + LSname + "_B5.TIF", nbrAr,
curFol + "Ind_NBR_" + LSname[:-5] + ".TIF")
nbrAr = []

if allIndices == 1 or allIndices == 8:
nbr2Ar = NBR2(reflectance6,reflectance7)
array_to_raster(curFol + LSname + "_B6.TIF", nbr2Ar,
curFol + "Ind_NBR2_" + LSname[:-5] + ".TIF")
nbr2Ar = []


#Save reflectance rasters to disk, arrays will be deleted either way
if saveReflectance:
for x in range(1,10):
exec( 'array_to_raster(curFol + LSname + "_B{0}.TIF", reflectance{0}, \
curFol + "Reflectance_B{0}.TIF")'.format(x) )
exec( 'reflectance{0} = []'.format(x) )
else:
for x in range(1,10):
exec( 'reflectance{0} = []'.format(x) )



#delete originally extracted files
if deleteOriginal == 2:
for x in orgList:
if x[-4] == ".":
os.remove(x)
else:
shutil.rmtree(x) #accounts for folder, os.remove only removes files
elif deleteOriginal == 3:
for x in orgList:
if x[-7:].lower() != "mtl.txt":
if x[-4] == ".":
os.remove(x)
else:
shutil.rmtree(x) #accounts for folder, os.remove only removes files






elif LSname[2] == '7' or LSname[2] == '5':
#Find the metafile and extract variables
for x in os.listdir(curFol):
#Following will only consider Landsat8 as valid input
if x[-7:].lower() == "mtl.txt":
#metaFile = tar.extractfile(x) #metaFiel is now a fileObject
metaFile = open(curFol + x, "r")
for lines in metaFile.readlines():
#remove lines' whitespaces
strLine = lines.replace(" ", "")

if strLine[:13] == 'SPACECRAFT_ID':
spacecraft_ID = strLine[15:-2]
if strLine[:8] == 'WRS_PATH':
wrs_path = strLine[9:]
if strLine[:7] == 'WRS_ROW':
wrs_row = strLine[8:]
if strLine[:11] == 'CLOUD_COVER':
cloud_cover = float(strLine[12:-1])
if strLine[:13] == 'DATE_ACQUIRED':
date_aquired_str = strLine[14:-1]
if strLine[:17] == 'SCENE_CENTER_TIME':
scene_center_time_str = strLine[18:-1]
if strLine[:8] == 'UTM_ZONE':
utm_zone = strLine[9:-1]

if strLine[:11] == 'SUN_AZIMUTH':
sun_azimuth = float(strLine[12:])
if strLine[:13] == 'SUN_ELEVATION':
sun_elevation = float(strLine[14:])

if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_1':
radiance_max_B1 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_1':
radiance_min_B1 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_2':
radiance_max_B2 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_2':
radiance_min_B2 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_3':
radiance_max_B3 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_3':
radiance_min_B3 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_4':
radiance_max_B4 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_4':
radiance_min_B4 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_5':
radiance_max_B5 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_5':
radiance_min_B5 = float(strLine[24:-2])
if LSname[2] == '7':
if strLine[:30] == 'RADIANCE_MAXIMUM_BAND_6_VCID_1':
radiance_max_B61 = float(strLine[31:-2])
if strLine[:30] == 'RADIANCE_MINIMUM_BAND_6_VCID_1':
radiance_min_B61 = float(strLine[31:-2])
if strLine[:30] == 'RADIANCE_MAXIMUM_BAND_6_VCID_2':
radiance_max_B62 = float(strLine[31:-2])
if strLine[:30] == 'RADIANCE_MINIMUM_BAND_6_VCID_2':
radiance_min_B62 = float(strLine[31:-2])
if LSname[2] == '5':
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_6':
radiance_max_B6 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_6':
radiance_min_B6 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_7':
radiance_max_B7 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_7':
radiance_min_B7 = float(strLine[24:-2])
if LSname[2] == '7':
if strLine[:23] == 'RADIANCE_MAXIMUM_BAND_8':
radiance_max_B8 = float(strLine[24:-2])
if strLine[:23] == 'RADIANCE_MINIMUM_BAND_8':
radiance_min_B8 = float(strLine[24:-2])

if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_1':
quant_max_B1 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_1':
quant_min_B1 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_2':
quant_max_B2 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_2':
quant_min_B2 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_3':
quant_max_B3 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_3' :
quant_min_B3 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_4' :
quant_max_B4 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_4' :
quant_min_B4 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_5':
quant_max_B5 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_5' :
quant_min_B5 = int(strLine[24:])
if LSname[2] == '7':
if strLine[:30] == 'QUANTIZE_CAL_MAX_BAND_6_VCID_1':
quant_max_B61 = int(strLine[31:])
if strLine[:30] == 'QUANTIZE_CAL_MIN_BAND_6_VCID_1':
quant_min_B61 = int(strLine[31:])
if strLine[:30] == 'QUANTIZE_CAL_MAX_BAND_6_VCID_2':
quant_max_B62 = int(strLine[31:])
if strLine[:30] == 'QUANTIZE_CAL_MIN_BAND_6_VCID_2':
quant_min_B62 = int(strLine[31:])
if LSname[2] == '5':
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_6':
quant_max_B6 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_6':
quant_min_B6 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_7':
quant_max_B7 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_7':
quant_min_B7 = int(strLine[24:])
if LSname[2] == '7':
if strLine[:23] == 'QUANTIZE_CAL_MAX_BAND_8' :
quant_max_B8 = int(strLine[24:])
if strLine[:23] == 'QUANTIZE_CAL_MIN_BAND_8' :
quant_min_B8 = int(strLine[24:])

if strLine[:20] == 'RADIANCE_MULT_BAND_1':
radiance_mult_B1 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_2':
radiance_mult_B2 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_3':
radiance_mult_B3 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_4':
radiance_mult_B4 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_5':
radiance_mult_B5 = float(strLine[21:])
if LSname[2] == '7':
if strLine[:27] == 'RADIANCE_MULT_BAND_6_VCID_1':
radiance_mult_B61 = float(strLine[28:])
if strLine[:27] == 'RADIANCE_MULT_BAND_6_VCID_2':
radiance_mult_B62 = float(strLine[28:])
if LSname[2] == '5':
if strLine[:20] == 'RADIANCE_MULT_BAND_6':
radiance_mult_B6 = float(strLine[21:])
if strLine[:20] == 'RADIANCE_MULT_BAND_7':
radiance_mult_B7 = float(strLine[21:])
if LSname[2] == '7':
if strLine[:20] == 'RADIANCE_MULT_BAND_8':
radiance_mult_B8 = float(strLine[21:])

if strLine[:19] == 'RADIANCE_ADD_BAND_1':
radiance_add_B1 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_2':
radiance_add_B2 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_3':
radiance_add_B3 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_4':
radiance_add_B4 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_5':
radiance_add_B5 = float(strLine[20:])
if LSname[2] == '7':
if strLine[:26] == 'RADIANCE_ADD_BAND_6_VCID_1':
radiance_add_B61 = float(strLine[27:])
if strLine[:26] == 'RADIANCE_ADD_BAND_6_VCID_2':
radiance_add_B62 = float(strLine[27:])
if LSname[2] == '5':
if strLine[:19] == 'RADIANCE_ADD_BAND_6':
radiance_add_B6 = float(strLine[20:])
if strLine[:19] == 'RADIANCE_ADD_BAND_7':
radiance_add_B7 = float(strLine[20:])
if LSname[2] == '7':
if strLine[:19] == 'RADIANCE_ADD_BAND_8':
radiance_add_B8 = float(strLine[20:])

metaFile.close()
members = []

#calculate radiation and toa relfectance
if LSname[2] == '7':
lastBand = 9
else:
lastBand = 8

for x in range(1,lastBand):
#read bands as arrays, bands 6 need to be handled separately
if x != 6:
exec( 'arrayDN{0} = singleTifToArray( \
curFol + LSname + "_B{0}.TIF")'.format(x) )

#convert to radiance and convert to 32-bit floating point for memory saving
exec( 'lambda{0} = ((radiance_max_B{0} - radiance_min_B{0}) / \
(quant_max_B{0} - quant_min_B{0})) * \
(arrayDN{0} - quant_min_B{0}) + \
radiance_min_B{0}'.format(x) )

exec( 'lambda{0} = lambda{0}.astype(np.float32)'.format(x))

#convert to reflectance and convert to 32-bit floating point for memory saving
esun = [1970,1842,1547,1044,225.7,82.06,1369][x-2] #band depending constant
e_s_dist = ES_dist(julDay)
exec( 'reflectance{0} = ( np.pi * lambda{0} * e_s_dist**2 ) / \
(esun * \
math.sin(math.radians(sun_elevation)))'.format(x) )

exec( 'reflectance{0} = reflectance{0}.astype(np.float32)'.format(x))

#empty first array for memory saving
exec( 'arrayDN{0} = []'.format(x) )


else:
if calcKelvin > 0:
if LSname[2] == '7':
if calcKelvin == 1:
startNum = 1
endNum = 3
elif calcKelvin == 2:
startNum = 1
endNum = 2
else:
startNum = 2
endNum = 3

for y in range(startNum,endNum):
exec( 'arrayDN{0}{1} = singleTifToArray( \
curFol + LSname + "_B{0}_VCID_{1}.TIF")'.format(x,y) )

exec( 'lambda{0}{1} = (( \
radiance_max_B{0}{1} - radiance_min_B{0}{1})/ \
(quant_max_B{0}{1} - quant_min_B{0}{1})) * \
(arrayDN{0}{1} - quant_min_B{0}{1}) + \
radiance_min_B{0}{1}'.format(x,y) )

exec( 'lambda{0}{1} = \
lambda{0}{1}.astype(np.float32)'.format(x,y))

k1_const_B61 = 666.09 #607.76 for LS5
k2_const_B62 = 1282.71 #1260.56 for LS5

exec ('t{0}{1} = k2_const_B62 / \
( np.log( (k1_const_B61 / lambda{0}{1}) + 1) )'.format(x,y) ) #T in Kelvin
exec ('t{0}{1} = t{0}{1}.astype(np.float32)'.format(x,y) )

exec( 'arrayDN{0}{1} = []'.format(x,y) )

exec( 'array_to_raster(curFol + LSname + \
"_B{0}_VCID_{1}.TIF", t{0}{1}, curFol + \
"Temp_B{0}{1}.TIF")'.format(x,y) )

exec( 't{0}{1} = []'.format(x,y) )


if LSname[2] == '5':
arrayDN6 = singleTifToArray(
curFol + LSname + "_B6.TIF")

lambda6 = ( ( (radiance_max_B6- radiance_min_B6)/
(quant_max_B6 - quant_min_B6) ) *
(arrayDN6 - quant_min_B6) +
radiance_min_B6 )
lambda6 = lambda6.astype(np.float32)

k1_const_B6 = 607.76
k2_const_B6 = 1260.56

t6 = k2_const_B6 / np.log( (k1_const_B6 / lambda6) + 1)
t6 = t6.astype(np.float32)
arrayDN6 = []
array_to_raster(curFol + LSname +
"_B6.TIF", t6, curFol + "Temp_B6.TIF")
t6 = []

continue


#Save radiation rasters to disk, arrays will be deleted either way
if saveRadiation:
for x in range(1,lastBand):
if LSname[2] == '7':
if x != 6:
exec( 'array_to_raster(curFol + LSname + "_B{0}.TIF", lambda{0}, \
curFol + "Radiation_B{0}.TIF")'.format(x) )
exec( 'lambda{0} = []'.format(x) )
else:
for y in range(1,3):
exec( 'array_to_raster(curFol + LSname + \
"_B{0}_VCID_{1}.TIF", lambda{0}{1}, curFol + \
"Radiation_B{0}{1}.TIF")'.format(x,y) )

exec( 'lambda{0}{1} = []'.format(x,y) )
if LSname[2] == '5':
exec( 'array_to_raster(curFol + LSname + "_B{0}.TIF", lambda{0}, \
curFol + "Radiation_B{0}.TIF")'.format(x) )
exec( 'lambda{0} = []'.format(x) )

else:
for x in range(1,lastBand):
if LSname[2] == '7':
if x != 6:
exec( 'lambda{0} = []'.format(x) )
else:
for y in range(1,3):
exec( 'lambda{0}{1} = []'.format(x,y) )

if LSname[2] == '5':
exec( 'lambda{0} = []'.format(x) )



#Calculate spectral indices
if allIndices > 0:
if allIndices == 1 or allIndices == 2:
ndviAr = NDVI(reflectance3,reflectance4)
array_to_raster(curFol + LSname + "_B4.TIF", ndviAr,
curFol + "Ind_NDVI_" + LSname[:-5] + ".TIF")
ndviAr = []

if allIndices == 1 or allIndices == 3:
eviAr = EVI(reflectance1,reflectance3,reflectance4,L=1)
array_to_raster(curFol + LSname + "_B3.TIF", eviAr,
curFol + "Ind_EVI_" + LSname[:-5] + ".TIF")
eviAr = []

if allIndices == 1 or allIndices == 4:
saviAr = SAVI(reflectance3,reflectance4,L=0.5)
array_to_raster(curFol + LSname + "_B4.TIF", saviAr,
curFol + "Ind_SAVI_" + LSname[:-5] + ".TIF")
saviAr = []

if allIndices == 1 or allIndices == 5:
msaviAr = MSAVI(reflectance3,reflectance4)
array_to_raster(curFol + LSname + "_B4.TIF", msaviAr,
curFol + "Ind_MSAVI_" + LSname[:-5] + ".TIF")
msaviAr = []

if allIndices == 1 or allIndices == 6:
ndmiAr = NDMI(reflectance4,reflectance5)
array_to_raster(curFol + LSname + "_B5.TIF", ndmiAr,
curFol + "Ind_NDMI_" + LSname[:-5] + ".TIF")
ndmiAr = []

if allIndices == 1 or allIndices == 7:
nbrAr = NBR(reflectance4,reflectance7)
array_to_raster(curFol + LSname + "_B5.TIF", nbrAr,
curFol + "Ind_NBR_" + LSname[:-5] + ".TIF")
nbrAr = []

if allIndices == 1 or allIndices == 8:
nbr2Ar = NBR2(reflectance5,reflectance7)
array_to_raster(curFol + LSname + "_B5.TIF", nbr2Ar,
curFol + "Ind_NBR2_" + LSname[:-5] + ".TIF")
nbr2Ar = []

#Save reflectance rasters to disk, arrays will be deleted either way
if saveReflectance:
for x in range(1,lastBand):
if x != 6:
exec( 'array_to_raster(curFol + LSname + "_B{0}.TIF", \
reflectance{0}, curFol + "Reflectance_B{0}.TIF")'.format(x) )
exec( 'reflectance{0} = []'.format(x) )

else:
for x in range(1,lastBand):
if x != 6:
exec( 'reflectance{0} = []'.format(x) )


#delete originally extracted files
if deleteOriginal == 2:
for x in orgList:
if x[-4] == ".":
os.remove(x)
else:
shutil.rmtree(x) #accounts for folder, os.remove only removes files
elif deleteOriginal == 3:
for x in orgList:
if x[-7:].lower() != "mtl.txt":
if x[-4] == ".":
os.remove(x)
else:
shutil.rmtree(x) #accounts for folder, os.remove only removes files




# Actual Working Example

inFol = "...\myInFol\" # Folder containing several .tar.gz rasters
outFol = "...\myOutFol\" # Folder that should contain the outputs

# Append all Landsat archives to a list
fileList = []
for files in os.listdir(inFol):
fileList.append(inFol + files)

# Execute function for each Landsat image, regardless of Landsat No.
# Function arguments are defined at the top, above all function definitions
for inFilex in fileList:
Landsat(inFilex,outFol,deleteOriginal,saveRadiation,
saveReflectance,calcKelvin,allIndices)