# PURPOSE: Automate the computation of NDVI values from Landsat images
# retrieved from http://earthexplorer.usgs.gov.
# Radiance and reflectance formulas provided by the Landsat Handbook
# (http://landsathandbook.gsfc.nasa.gov/data_prod/prog_sect11_3.html)
#
library(raster)
### Function computes fraction of mean distance to sun (adopted from 'Landsat' R package)
Sdist = function (adate)
{
edist <- julian(as.Date(adate), origin = as.Date(paste(substring(adate,
1, 4), "12", "31", sep = "-")))[[1]]
edist <- 1 - 0.016729 * cos((2 * pi) * (0.9856 * (edist -
4)/360))
edist
}
### End function
# Some constants
E.red = 1551 # Red solar spectral irradiance
E.nir = 1044 # NIR solar spectral irradiance
# Get list of landsat metadata files in directory
flist <- list.files(pattern = "MTL\\.txt$")
# Loop through each file in list
for (i in 1:length(flist) ){
# Open metadata file
inFile = read.delim(flist[i],header=F,sep="=", stringsAsFactors=F)
inFile$V1 = gsub(' +', '', inFile$V1) # Remove spaces from 1st column
# Extract parameters from metadata file
sun.elev = as.numeric(inFile$V2[inFile$V1 == "SUN_ELEVATION"]) * pi / 180 # Solar elevation in radians
lmax.red = as.numeric(inFile$V2[inFile$V1 == "LMAX_BAND3"])
lmin.red = as.numeric(inFile$V2[inFile$V1 == "LMIN_BAND3"])
lmax.nir = as.numeric(inFile$V2[inFile$V1 == "LMAX_BAND4"])
lmin.nir = as.numeric(inFile$V2[inFile$V1 == "LMIN_BAND4"])
Qmax.red = as.numeric(inFile$V2[inFile$V1 == "QCALMAX_BAND3"])
Qmin.red = as.numeric(inFile$V2[inFile$V1 == "QCALMIN_BAND3"])
Qmax.nir = as.numeric(inFile$V2[inFile$V1 == "QCALMAX_BAND4"])
Qmin.nir = as.numeric(inFile$V2[inFile$V1 == "QCALMIN_BAND4"])
L.date = inFile$V2[inFile$V1 == "ACQUISITION_DATE"]
L.date = gsub(' +', '', L.date) # Remove leading space
file.red = inFile$V2[inFile$V1 == "BAND3_FILE_NAME"]
file.nir = inFile$V2[inFile$V1 == "BAND4_FILE_NAME"]
# Open raster layers
red = raster(file.red)
nir = raster(file.nir)
# Compute Radiance (L)
red.L = (lmax.red - lmin.red)/(Qmax.red - Qmin.red) * (red - Qmin.red) + lmin.red
nir.L = (lmax.nir - lmin.nir)/(Qmax.nir - Qmin.nir) * (nir - Qmin.nir) + lmin.nir
# Compute Reflectance (R)
red.R = (pi * red.L * Sdist(L.date)^2) / (E.red * cos(sun.elev))
nir.R = (pi * nir.L * Sdist(L.date)^2) / (E.nir * cos(sun.elev))
# Compute NDVI
NDVI = (nir.R - red.R) / (nir.R + red.R)
# Save NDVI raster to file
outFile = paste("NDVI_",L.date,".img",sep="")
writeRaster(NDVI,filename=outFile,format="HFA", overwrite=F)
}