install.packages("rmarkdown")
install.packages("tidyverse", dependencies = TRUE)
install.packages("raster", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)

library(raster)

#blue
b2 <- raster('data/B2.tif')
#green
b3 <- raster('data/B3.tif')
#red
b4 <- raster('data/B4.tif')

landsatRGB <- stack(b4, b3, b2) #order is impt

plotRGB(landsatRGB, 
               stretch = "lin")

filenames <- paste0('data/B', 1:5, ".tif")

landsat <- stack(filenames)

names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')

crs(landsat)

library(sf)

sgshp <- st_read("data/boundary/sg.shp")

crs(sgshp)

sgshp <- st_transform(sgshp, crs = crs(landsat))

landsat <- crop(landsat, sgshp) # crop to rectangle

landsat <- mask(landsat, sgshp)

landsatRGB <- subset(landsat, c(4,3,2))
plotRGB(landsatRGB,
	stretch = "lin")

ndvi <- function(x, y) {
(x - y) / (x + y)
}

landsatNDVI <- overlay(landsat[[5]], landsat[[4]], fun = ndvi)

landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1
landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 1

plot(landsatNDVI, 
col = rev(terrain.colors(10)), 
main = "Landsat 8 NDVI",
axes = FALSE, box = FALSE)

hist(landsatNDVI,
main = "Distribution of NDVI values", xlab = "NDVI", 
xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)

landsatGreen <- reclassify(landsatNDVI, c(-1, 0.2, NA)) #-1 to 0.2 becomes NA
plot(landsatGreen, col = "darkgreen", 
axes = FALSE, box = FALSE, legend = FALSE)

landsatBlue <- reclassify(landsatNDVI, c(0, 1, NA)) # 0 to 1 becomes NA
plot(landsatBlue, col = "blue", 
axes = FALSE, box = FALSE, legend = FALSE,) 

reclass_m <- matrix(c(-Inf, 0, 1, #water
0, 0.2, 2, #urban
0.2, Inf, 3), #veg
ncol = 3, byrow = TRUE)

reclass_m

landsatCover <- reclassify(landsatNDVI, rcl = reclass_m)

plot(landsatCover,
col = c("blue", "grey", "darkgreen"),
legend = FALSE,
axes = FALSE,
box = FALSE,
main = "Land cover (mosaic) in Singapore")
legend("bottomright",
legend = c("Water", "Urban", "Vegetation"),
fill = c("blue", "grey", "darkgreen"),
border = FALSE,
bty = "n")


writeRaster(landsatCover, 
filename = "data/landsat_landcover.tif", 
overwrite = TRUE)


landsatCover <- raster('data/landsat_landcover.tif')

library(landscapemetrics)
library(landscapetools)

show_landscape(landsatCover, discrete = TRUE)

check_landscape(landsatCover)

list_lsm()

list_lsm(level = "patch", type = "area and edge metric")

lsm_p_area(landsatCover)

lsm_c_ca(landsatCover)

lsm_c_area_mn(landsatCover) #mean area

lsm_c_pland(landsatCover) #percentage of landscape

lsm_l_ta(landsatCover) #total area of landscape

calculate_lsm(landsatCover, 
what = c("lsm_l_ta", #total area of landscape
"lsm_c_ca", #total area of each class
"lsm_c_pland"), #proportional area of each class
full_name = TRUE) #include info about metrics in results

calculate_lsm(landsatCover, 
level = "class",
type = "aggregation metric",
full_name = TRUE) #include info about metrics in results

calculate_lsm(landsatCover, 
              level = "landscape",
              type = "diversity metric",
              full_name = TRUE)