library(raster)
library(tidyverse)
library(getlandsat)
library(sf)
library(mapview)
library(osmdata)
library(stats)

Question 1

uscities = read_csv("../data/uscities.csv")

bb = uscities %>%
  filter(city == "Palo") %>%
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_as_sf()

# mapview(bb)

Question 2

bbwgs = bb %>% st_transform(4326)
bb = st_bbox(bbwgs)
meta = read_csv("../data/palo-flood-scene.csv")

files = lsat_scene_files(meta$download_url) %>%
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file)) %>%
  arrange(file) %>%
  pull(file)

st = sapply(files, lsat_image)

s = stack(st) %>%
  setNames(c(paste0("band", 1:6)))

cropper = bbwgs %>% st_transform(crs(s))

r = crop(s, cropper)
dimensions : 340, 346, 117640, 6 (nrow, ncol, ncell, nlayers)
crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
resolution : 30, 30 (x, y)

Question 3

par(mfrow = c(2, 2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
plotRGB(r, r = 5, g = 4, b = 3, stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4, stretch = "hist")
plotRGB(r, r = 6, g = 3, b = 2, stretch = "lin")

The color stretch adds contrast to the plot to make the colors more enhanced.

Question 4

Step 1

ndvi = (r$band5 - r$band4) / (r$band5 + r$band4)
ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)
mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)
wri = (r$band3 + r$band4) / (r$band5 + r$band6)
swi  = (1 / (sqrt(r$band2 - r$band6)))

stack = stack(ndvi, ndwi, mndwi, wri, swi) %>%
  setNames(c("Normalized difference vegetation index",
             "Normalized Difference water Index",
             "Modified normalized difference water index",
             "Water ratio index",
             "Simple water index"))

palette = colorRampPalette(c("blue", "white", "red"))
plot(stack, col = palette(256))

Step 2

thresholding_ndvi = function(x){ifelse(x <= 0, 1, 0)}
thresholding_ndwi = function(x){ifelse(x >= 0, 1, 0)}
thresholding_mndwi = function(x){ifelse(x >= 0, 1, 0)}
thresholding_wri = function(x){ifelse(x >= 1, 1, 0)}
thresholding_swi = function(x){ifelse(x <= 5, 1, 0)}

flood_ndvi = calc(ndvi, thresholding_ndvi)
flood_ndwi = calc(ndwi, thresholding_ndwi)
flood_mndwi = calc(mndwi, thresholding_mndwi)
flood_wri = calc(wri, thresholding_wri)
flood_swi = calc(swi, thresholding_swi)

flood_stack = stack(flood_ndvi, flood_ndwi, flood_mndwi, flood_wri, flood_swi) %>%
  setNames(c("NDVI",
             "NDWI",
             "MNDWI",
             "WRI",
             "SWI"))

plot(flood_stack, col = c("white", "blue"))

Question 5

Step 1

set.seed(1)

Step 2

values = getValues(r)
dim(values)
## [1] 117640      6
# 117640    6

# There are 117,640 rows and 6 columns, meaning there is one column for each band of the raster stack

values = na.omit(values)

k12 = kmeans(values, centers = 12)
k6 = kmeans(values, centers = 6)
k3 = kmeans(values, centers = 3)

kmeans_raster12 = flood_stack$NDVI
kmeans_raster6 = flood_stack$NDVI
kmeans_raster3 = flood_stack$NDVI

values(kmeans_raster12) = k12$cluster
plot(kmeans_raster12)

values(kmeans_raster6) = k6$cluster
plot(kmeans_raster6)

values(kmeans_raster3) = k3$cluster
plot(kmeans_raster3)

Step 3

t = table(getValues(flood_swi), getValues(kmeans_raster12))

index = which.max(t[2,])
# 9

thresholding_index = function(x){ifelse(x == 9, 1, 0)}

flood_index = calc(kmeans_raster12, thresholding_index)

plot(flood_index)

flood_stack = flood_stack %>% addLayer(flood_index) %>% 
  setNames(c("NDVI",
             "NDWI",
             "MNDWI",
             "WRI",
             "SWI",
             "Index"))
plot(flood_stack, col = c("white", "blue"))

Question 6

rsum_plot = sum(flood_stack)

stats = cellStats(flood_stack, sum)
stats = (stats*900)/1000
knitr::kable(stats, caption = c("Area of Cells"), col.names = "Area")
Area of Cells
Area
NDVI 5999.4
NDWI 6490.8
MNDWI 10745.1
WRI 7622.1
SWI 13680.9
Index 7261.2 
plot(rsum_plot, col = blues9)

uncertainty = calc(flood_stack, fun = sum)
plot(uncertainty, col = blues9)

Some of the cell values aren’t even because the data is averages, meaning there is overlap in some of the cells.