## ----setup, echo=TRUE, include=FALSE------------------------------------------
library(knitr)
library(raster)


## ----lsstack------------------------------------------------------------------
library(raster)
raslist <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- stack(raslist)
landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]


## ----rs3vi--------------------------------------------------------------------
vi <- function(img, k, i) {
  bk <- img[[k]] 
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}


## ----rs3ndvi------------------------------------------------------------------
# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")


## ----rs3ndvialt---------------------------------------------------------------
vi2 <- function(x, y) {
	(x - y) / (x + y)
}
ndvi2 <- overlay(landsat[[5]], landsat[[4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")


## ----rs2hist------------------------------------------------------------------
# view histogram of data
hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))


## ----rs3veg1------------------------------------------------------------------
veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetation')


## ----rs3land------------------------------------------------------------------
land <- reclassify(ndvi, c(-Inf, 0.25, NA,  0.25, 0.3, 1,  0.3, Inf, NA))
plot(land, main = 'What is it?')


## ----rs3rgb1------------------------------------------------------------------
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")
plot(land, add=TRUE, legend=FALSE)


## ----rs3veg2------------------------------------------------------------------
vegc <- reclassify(ndvi, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI based thresholding')


## ----rs3pca1------------------------------------------------------------------
set.seed(1)
sr <- sampleRandom(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")


## ----rs3pca2------------------------------------------------------------------
pca <- prcomp(sr, scale = TRUE)
pca
screeplot(pca)

pci <- predict(landsat, pca, index = 1:2)
plot(pci[[1]])


## ----rs3rgb2, fig.width = 8, fig.height = 4-----------------------------------
pc2 <- reclassify(pci[[2]], c(-Inf,0,1,0,Inf,NA))
par(mfrow = c(1,2))
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plotRGB(landsatFCC, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin", main = "Landsat False Color Composite")
plot(pc2, legend = FALSE, add = TRUE)