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

#Question 1

bb <- read_csv("~/github/geog-176A-labs/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()
plot(bb)

#Question 2

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

meta = read_csv("~/github/geog-176A-labs/data/palo-flood.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(paste0("band", 1:6))
    plot(s)

s
## class      : RasterStack 
## dimensions : 7811, 7681, 59996291, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 518085, 748515, 4506585, 4740915  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 
## names      : band1, band2, band3, band4, band5, band6 
## min values :     0,     0,     0,     0,     0,     0 
## max values : 65535, 65535, 65535, 65535, 65535, 65535

dimensions : 7811, 7681, 59996291, 6 (nrow, ncol, ncell, nlayers)

resolution : 30, 30 (x, y)

crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

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

r = crop(s, cropper)
r = r %>% 
  setNames(c("Coastal Aerosol", "Blue", "Green", "Red", "Near Infrared", "SWIR 1"))

plot(r)

r
## class      : RasterBrick 
## dimensions : 340, 346, 117640, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594075, 604455, 4652535, 4662735  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : Coastal.Aerosol,  Blue, Green,   Red, Near.Infrared, SWIR.1 
## min values :            8566,  7721,  6670,  6057,          5141,   4966 
## max values :           18482, 19079, 19592, 21965,         24772,  23896

dimensions : 340, 346, 117640, 6 (nrow, ncol, ncell, nlayers)

resolution : 30, 30 (x, y)

crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

#Question 3

R-G-B (Natural Color Image)

plotRGB(r, r = 4, g = 3, b = 2)

NIR-R-G(fa) (Traditional Color Infrared (CIR) Image)

plotRGB(r, r = 5, g = 4, b = 3)

NIR-SWIR1-R (False Color good for picking out land from water)

plotRGB(r, r = 5, g = 6, b = 4)

SWIR2-SWIR1-R (False Color useful for visualizing urban environments)

plotRGB(r, r = 7, g = 6, b = 4)

Colored stretch is an effective and simple way to keep a visual track of different products because colored stretch is easier for people to quickly identify and separate different elements, ensuring a specific element isn’t mixed up with others.

RGB image comparison w stretch=“hist”

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

RGB image comparison w stretch=“lin”

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 = "lin")
plotRGB(r, r = 7, g = 6, b = 4, stretch = "lin")

#Question 4

ndvi = (r$Near.Infrared - r$Red) / (r$Near.Infrared + r$Red)
ndwi = (r$Green - r$Near.Infrared) / (r$Green + r$Near.Infrared)
mndwi = (r$Green - r$SWIR.1) / (r$Green + r$SWIR.1)
wri = (r$Green + r$Red) / (r$Near.Infrared + r$SWIR.1)
swi = 1 / (sqrt(r$Blue - r$SWIR.1))
stack = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue","white","red"))
plot(stack, col = palette(256))

All 5 images emphasize the water feature, and all of them are in blue and red. However, the water features are presented in different color. In NDVI and SWI, the threshold of water features are less than 0. In NDWI, MNDWI, and WIR, the threshold of water features are greater than 0. The range of the water threshold is the greatest in MNDWI and the least in SWI

ndvi2 = function(x){ifelse(x <= 0,1, 0)}
ndwi2 = function(x){ifelse(x >= 0,1, 0)}
mndwi2 = function(x){ifelse(x >= 0,1, 0)}
wri2 = function(x){ifelse(x <= 1,1, 0)}
swi2 = function(x){ifelse(x<=5, 1, 0)}

ndvi3 = calc(ndvi, ndvi2)
ndwi3 = calc(ndwi, ndwi2)
mndwi3 = calc(mndwi, mndwi2)
wri3 = calc(wri, wri2)
swi3 = calc(swi, swi2)

threshold = stack(ndvi3, ndwi3, mndwi3, wri3, swi3) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI")) 
palette = colorRampPalette(c("white","blue"))
plot(threshold, col = palette(256))

#Question 5

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

There are 117,640 rows, 6 columns.

k12 = kmeans(values, 12, iter.max = 100)
kmeans_raster = r$Coastal.Aerosol
values(kmeans_raster) = k12$cluster
plot(kmeans_raster)

kmeans_vals = getValues(kmeans_raster)
binary_vals = getValues(ndvi3)
table <- table(binary_vals,kmeans_vals)

which.max(table)
## [1] 13
t_kmeans = function(x){ifelse(x == 3, 1, 0)}
f_kmeans = calc(kmeans_raster, t_kmeans)

threshold = addLayer(threshold, f_kmeans) %>% 
  setNames(c("ndvi", "ndwi", "mndwi", "wri", "swi", "kmeans"))
palette = colorRampPalette(c("white","blue"))
plot(threshold, col = palette(256))

#Question 6

kabletable = cellStats(threshold, sum)
knitr::kable(kabletable, caption = "Number of Flooded Cells", col.names = c("Number"))
Number of Flooded Cells
Number
ndvi 6666
ndwi 7212
mndwi 11939
wri 109172
swi 15201
kmeans 7395 
areakable = kabletable * 900
knitr::kable(areakable, caption = "Area of Flooded Cells (m^2)", col.names = c("Area"))
Area of Flooded Cells (m^2)
Area
ndvi 5999400
ndwi 6490800
mndwi 10745100
wri 98254800
swi 13680900
kmeans 6655500

Extra Creedit

point = st_point(c(-91.78948,42.06306)) %>% 
  st_sfc(crs = 4326) %>% 
  st_transform(crs(threshold)) %>%
  as_Spatial()
print("-91.78948,42.06306")
## [1] "-91.78948,42.06306"
raster::extract(threshold, point)
##      ndvi ndwi mndwi wri swi kmeans
## [1,]    1    1     1   0   1      0