R: Compute NDVI from Landsat images

18 February, 2012 (23:44) | R, Remote Sensing, Spatial analysis | By: Manuel Gimond

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