Libraries

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

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()

mapview(bb)

Question 2

bbwgs= bb %>% st_transform(4326)

bb=st_bbox(bbwgs)


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

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

Reading in created .csv data

md = read_csv("../data/palo_flood.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 Color

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")

Color Infrared

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

False Color Water Focus

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

My Choice

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

Stretching determines how the values are represented on the image. Histogram stretching tends to be lighter, while linear stretching seems to be more discernible.

Question 4

Creating the 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 masks

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)

Creating and plotting 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)

mapview(flood2)

Question 5

Creating the Kmeans cluster

set.seed(1)

#1. Splited Raster and get data to cluster
r_extract = getValues(r)
r_nna = na.omit(r_extract) #No NAs, no need to omit
#data was extracted band by band (there are 6 columns)

#2. Calcualting the KMeans Cluster
kmeans_r = kmeans(r_extract, 12)
fviz_cluster(kmeans_r, geom="point", data = r_extract)

#3. Applying the cluster to raster structure
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?

After finding the lat and long in the video:

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 our area as a flood, meaning some of the rasters missed marking this area as a flood.