Below I have gathered Landsat data which will be used to calculate NVDI:

all_landsat_bands <- list.files("C:/Users/kaitlyn/documents/earth-analytics/data/week6/Landsat/LC80340322016205-SC20170127160728/crop",
           pattern=glob2rx("*band*.tif$"),
           full.names = T)
landsat_stack_csf <- stack(all_landsat_bands)

The script above instructs r to gather all the bands and stack them together.

The red and infrared bands will be used to calculate the index using a specified ratio:

landsat_ndvi <- (landsat_stack_csf[[5]] - landsat_stack_csf[[4]]) / (landsat_stack_csf[[5]] + landsat_stack_csf[[4]])

Then the result is plotted:

plot(landsat_ndvi,
     main="Landsat derived NDVI\n 23 July 2016")

The values are examined in the histogram below: