Libraries

# SPDS
library(tidyverse)
library(sf)
library(units)

# Data
library(USAboundaries)
library(rnaturalearth)

# Visualization
library(gghighlight)
library(ggrepel)
library(knitr)
library(kableExtra)
library(readxl)
library(leaflet)
library(raster) # Raster Data handling
library(tidyverse) # Data Manipulation
library(getlandsat) # keyless Landsat data (2013-2017)
library(sf) # Vector data processing
library(mapview) # Rapid Interactive visualization
library(factoextra)

Question 1: AOI Identification

Palo = 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

Step 2: Read in Metadata

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)

Step 3:

st = sapply(files, lsat_image)

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

The dimensions of the stacked image are 7811 rows by 7681 columns. The Coordinate Reference System used is UTM Zone 15 and the cell resolution is 30 x 30.

Step 4

crop = Palo %>%
  st_transform(crs(s))

The dimensions of the stacked image are 340 rows by 346 columns. The Coordinate Reference System used is UTM Zone 15 and the cell resolution is 30 x 30.

Question 3 – RGB Images

r = crop(s, crop)
par(mfrow = c(2,2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
#natural
plotRGB(r, r = 5, g = 4, b = 3, stretch = "lin")
#Traditional Color Infrared
plotRGB(r, r = 5, g = 6, b = 4, stretch = "lin")
#False Color used for Water focus
plotRGB(r, r = 5, g = 7, b = 1, stretch = "lin")

#uses NIR, SWIR2, Coastal Aerosols used to visualize vegetation in water

Part 2

Applying a color stretch provides greater contrast between the coloration of cells within the raster, and can make it easier to observe desired spatial trends.

Question 4

palette = colorRampPalette(c("blue", "white", "red"))
ndvi = (r$band5 - r$band4) / (r$band5 + r$band4)
plot(ndvi, col = palette(256))

ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)
plot(ndwi, col = palette(256))

mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)
plot(mndwi, col = palette(256))

wri = (r$band3 + r$band4) / (r$band5 + r$band6)
plot(wri, col = palette(256))

swi = (1 / (sqrt(r$band2 - r$band6)))
plot(swi, col = palette(256))

index_stack = stack(ndvi, ndwi, mndwi, wri, swi) %>%
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(index_stack, col = palette(256))

Each of the plots help visualize slightly different components of the same area of interest. The NDVI plot shows the vegetated land in red, and non-vegetated land in blue. The NDWI plot shows is a similar plot, but instead shows areas which are water in red and dry-land areas in blue. The MNDWI plot enhances the visualization of the NDWI plot, showing further distinction between water (red) and land (blue). The WRI plot goes even further, normalizing the colors of water and land regions. The final plot completely removes land (shows in white) and leaves us to visualize only the water areas in blue.

Step 2 Thresholding

thresholding1 = function(x){ifelse(x <= 0, 1, NA)}
flood1 = calc(ndvi, thresholding1)
thresholding2 = function(x){ifelse(x >= 0, 1, NA)}
flood2 = calc(ndwi, thresholding2)
thresholding3 = function(x){ifelse(x >= 0, 1, NA)}
flood3 = calc(mndwi, thresholding3)
thresholding4 = function(x){ifelse(x >= 1, 1, NA)}
flood4 = calc(wri, thresholding4)
thresholding5 = function(x){ifelse(x <= 5, 1, NA)}
flood5 = calc(swi, thresholding5)
flood_stack = stack(flood1, flood2, flood3, flood4, flood5) %>%
  setNames(c("NDVI Flood", "NDWI Flood", "MNDWI Flood", "WRI Flood", "SWI Flood"))
plot(flood_stack, col = "blue")

Question 5

set.seed(09102020)
getValues(r) %>%
  dim() %>%
  na.omit()
## [1] 117640      6

The number of extracted values matches the number of cells in the raster, or the dimensions multiplied together. This means that the values were extraceted cell by cell, and that each cell was assigned a single value.

data_r = getValues(r) %>%
  na.omit()
kmean = kmeans(data_r, 12)
fviz_cluster(kmean, geom="point", data = data_r)

thresholding = function(x){ifelse(x <= 0, 1, 0)}
findflood = calc(ndvi, thresholding)
kmeansraster = findflood 
values(kmeansraster) = kmean$cluster
plot(kmeansraster, col = viridis::viridis(12))

com_table = table(values(findflood),values(kmeansraster))
kmeansraster[kmeansraster != which.max(com_table[2,])] = 0
kmeansraster[kmeansraster != 0] = 1
plot(kmeansraster)

thresholdings2 = function(x){ifelse(x >= 0, 1, 0)}
findflood2 = calc(ndwi, thresholdings2)
thresholdings3 = function(x){ifelse(x >= 0, 1, 0)}
findflood3 = calc(mndwi, thresholdings3)
thresholdings4 = function(x){ifelse(x >= 1, 1, 0)}
findflood4 = calc(wri, thresholdings4)
thresholdings5 = function(x){ifelse(x <= 5, 1, 0)}
findflood5 = calc(swi, thresholdings5)
stackfindflood = stack(findflood, findflood2, findflood3, findflood4, findflood5, kmeansraster)
plot(stackfindflood)

Question 6

sumsff = sum(stackfindflood)
plot(sumsff, col = blues9)

(cellStats(stackfindflood, sum) * res(sumsff)^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.5780

Extra Credit

AOI = st_point(c(-91.78967, 42.06290)) %>%
  st_sfc(crs = 4326) %>%
  st_as_sf() %>%
  st_transform(st_crs(stackfindflood)) 
raster::extract(sumsff, AOI)
## [1] 2

Two of the rasters identified the point as a flooding area.