## ----lsstack------------------------------------------------------------------
library(terra)
rfiles <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(rfiles)
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, fig.height=4.23-------------------------------------------------
# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")


## ----rs3ndvialtfun------------------------------------------------------------
vi2 <- function(x, y) {
	(x - y) / (x + y)
}


## ----rs3ndvialt, fig.height=4.23----------------------------------------------
nir <- landsat[[5]]
red <- landsat[[4]]
ndvi2 <- lapp(c(nir, red), fun = vi2)

# or in one line
#ndvi2 <- lapp(landsat[[5:4]], fun=vi2)

plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")


## ----rs2hist------------------------------------------------------------------
hist(ndvi, main = "NDVI values", xlab = "NDVI", ylab= "Frequency",  
   	col = "wheat", xlim = c(-0.5, 1),  breaks = 30, xaxt = "n")

axis(side=1, at = seq(-0.6, 1, 0.2), labels = seq(-0.6, 1, 0.2))



## ----rs3veg1, fig.height=4.23-------------------------------------------------
veg <- clamp(ndvi, 0.4, values=FALSE)
plot(veg, main='Vegetation')


## ----rs3land, fig.height=4.23-------------------------------------------------
m <- c(-Inf, 0.25, NA,  0.25, 0.3, 1,  0.3, Inf, NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
land <- classify(ndvi, rcl)
plot(land, main = 'What is it?')


## ----rs3rgb1, fig.height=4.23-------------------------------------------------
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin")
plot(land, add=TRUE, legend=FALSE)


## ----rs3veg2, fig.height=4.23-------------------------------------------------
m <- c(-1,0.25, 0.3, 0.4, 0.5, 1)
vegc <- classify(ndvi, m)
plot(vegc, col = rev(terrain.colors(4)), main = 'NDVI based thresholding')


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


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


## ----rs3pca2b, fig.width=8, fig.height=3--------------------------------------
pca_predict2 <- function(model, data, ...) {
  predict(model, data, ...)[,1:2]
}
pci <- predict(landsat, pca, fun=pca_predict2)
plot(pci)


## ----rs3rgb2------------------------------------------------------------------
# quick look at the histogram of second component
hist <- pci[[2]]
m <- c(-Inf,-3,NA, -3,-2,0, -2,-1,1, -1,0,2, 0,1,3, 1,2,4, 2,6,5, 6,Inf,NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
rcl
pcClass <- classify(pci[[2]], rcl)


## ----rs3rgb2plot, fig.width=10, fig.height=4----------------------------------
par(mfrow=c(1,2))
plotRGB(landsatFCC, stretch = "lin", main="False Color", mar=c(3.1, 3.1, 2.1, 2.1))
plot(pcClass, main="PCA")