library(raster)
library(tidyverse)
library(getlandsat) 
library(sf) 
library(mapview)
library(osmdata)
library(factoextra)
library(leaflet)

#Question 1

bb = read_csv('../data/uscities.csv') %>% 
  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()

#Question 2

bbwgs = bb %>% st_transform(4326)
bb = st_bbox((bbwgs))

q2_image = getlandsat::lsat_scenes() %>% 
  filter(min_lat <= bb$ymin, min_lon <= bb$xmin) %>% 
  filter(max_lat >= bb$ymax, max_lon >= bb$xmax) %>% 
  filter(as.Date(acquisitionDate) == as.Date("2016-09-26"))

write.csv(q2_image, file = "../data/palo-flood-scene.csv")

Step 2

md = read_csv("../data/palo-flood-scene.csv")

files = lsat_scene_files(md$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)

#Question 3 True

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

Infrared

plotRGB(r, r=5 ,g=4 ,b=3, stretch= "lin")

plotRGB(r, r=5 ,g=4 ,b=3, stretch= "hist")

False Color

plotRGB(r, r=5 ,g=6 ,b=4, stretch= "lin")

plotRGB(r, r=5 ,g=6 ,b=4, stretch= "hist")

My Choice

plotRGB(r, r=2 ,g=4 ,b=6, stretch= "lin")

plotRGB(r, r=2 ,g=4 ,b=6, stretch= "hist")

Stretching determines the representation of each value. Histogram stretching tends to be lighter and brighter, while linear stretching seems to be darker

Question 4

5 Different Masks

palette= colorRampPalette(c("blue","white","red"))

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)

Plotting

plot(ndvi,col=palette(256), legend=FALSE)

plot(ndwi,col=palette(256), legend=FALSE)

plot(mndwi,col=palette(256), legend=FALSE)

plot(wri,col=palette(256), legend=FALSE)

plot(swi,col=palette(256), legend=FALSE)

Water Thresholds

thresholding1= function(x){ifelse(x <= 0,1,NA)}
thresholding2= function(x){ifelse(x >= 0,1,NA)}
thresholding3= function(x){ifelse(x >= 0,1,NA)} 
thresholding4= function(x){ifelse(x >= 1,1,NA)}
thresholding5= function(x){ifelse(x <= 5,1,NA)}

flood1= calc(ndvi,thresholding1)
flood2= calc(ndwi,thresholding2)
flood3= calc(mndwi,thresholding3)
flood4= calc(wri,thresholding4)
flood5= calc(swi,thresholding5)

plot(flood1, col= "blue", legend=FALSE)

plot(flood2, col= "blue", legend=FALSE)

plot(flood3, col= "blue", legend=FALSE)

plot(flood4, col= "blue", legend=FALSE)

plot(flood5, col= "blue", legend=FALSE)

Question 5

Kmeans

set.seed(1)

r_extract = getValues(r)
r_nna = na.omit(r_extract) 


kmeans_r = kmeans(r_extract, 12)
fviz_cluster(kmeans_r, geom="point", data = r_extract)

thresholding_s= function(x){ifelse(x >= 0,1,0)}
flood_s= calc(ndwi,thresholding_s)
kmeans_raster = flood_s
values(kmeans_raster) = kmeans_r$cluster
plot(kmeans_raster, col = viridis::viridis(12))

Finding the correct KMeans cluster that represents our water pixels

com_table = table(values(flood_s), values(kmeans_raster))

which.max(com_table[2,])
## 12 
## 12
kmeans_raster[kmeans_raster != which.max(com_table[2,])] = 0
kmeans_raster[kmeans_raster != 0] =1
plot(kmeans_raster)

Changing remaining NAs to 0 for analysis

thresholdings1= function(x){ifelse(x <= 0,1,0)}
thresholdings3= function(x){ifelse(x >= 0,1,0)} 
thresholdings4= function(x){ifelse(x >= 1,1,0)}
thresholdings5= function(x){ifelse(x <= 5,1,0)}

floods1= calc(ndvi,thresholdings1)
floods3= calc(mndwi,thresholdings3)
floods4= calc(wri,thresholdings4)
floods5= calc(swi,thresholdings5)

fs = stack(floods1, flood_s, floods3, floods4, floods5, kmeans_raster)
plot(fs)

sum_fs = fs %>% sum() 
plot(sum_fs, col = blues9)

(cellStats(fs, sum) * res(fs)^2) /1e6
## layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 
##  5.9994  6.4908 10.7451  7.6221 13.6809  7.9371

Extra Credit

How many layers classified the marked area in the video as a flood?

aoi = st_point(c(-91.78967, 42.06290)) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf() %>% 
  st_transform(st_crs(fs))

raster::extract(sum_fs, aoi)
## [1] 3

There are 3 layers that classified this area as a flood; it is assumed that some rasters missed the marking of the area as a flood