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

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:

meta= read_csv("../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(c(paste0("band", 1:6)))

The dimensions of my stacked image are: 7811 (nrow)/ 7681 (ncol)/ 59996291 (ncell)/ 6 (nlayers)

The CRS is: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

The cell resolution is: 30, 30 (x, y)

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

r=crop(s, cropper)

The dimensions of my cropped image stacked are: 340 (nrow)/ 346 (ncol)/ 117640 (ncell)/ 6 (nlayers)

The CRS is: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

The cell resolution is: 30, 30 (x, y)

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 Infared

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

Color stretching determines how the values are represented on the image.

Question 4:

Step 1

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)

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)

All five images give us the results of different index methods which we first had to state. The normalized methods in short minimizes the effects of illumination and enhance spectral features that aren’t originally visible. Our water threshold which is given show us areas that both retain and don’t retain water & for NDVI index, numbers less than 0 = water shown in blue and numbers greater than 0 = vegetation shown in red. For both NDWI and MNDWI Index’s, numbers greater than 0 = water shown in red and numbers less than 0 = vegetation without water presence. For WRI Index, it shows the water ratio by displaying water presence in cells with numbers greater than 1 in red. The same can be said for the SWI Index since it displays water presence in cells with values greater than 5.

Step 2

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)

changing NA’s to 0s for Q5
thresholdings1= function(x){ifelse(x <= 0,1,0)}
thresholdings2= 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)
floods2= calc(ndwi,thresholdings2) 
floods3= calc(mndwi,thresholdings3)
floods4= calc(wri,thresholdings4)
floods5= calc(swi,thresholdings5)

Question 5:

Genaric Steps of KMeans Clustering

  1. Split raster data/structure into just data
  2. Calculate KMeans clustering on data
  3. Add clustered data to original raster structure
set.seed(1)

Must be set for any random process in R so that results are the same.

1. Spliting Raster and getting data to cluster

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

There were no NAs, no need to use na.omit() and 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

kmeans_raster = floods2
values(kmeans_raster) = kmeans_r$cluster
plot(kmeans_raster, col = viridis::plasma(12))

Step 5.3

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

which.max(com_table[2,])
## 12 
## 12

Here we determined which group in our KMeans clustering correlated to water pixels in our raster.

rast_thresh = function(x){ifelse(x != which.max(com_table[2,]),0,1)}
kmeans_floodmask = calc(kmeans_raster, rast_thresh)

floodst = stack(floods1, floods2, floods3, floods4, floods5, kmeans_floodmask)
plot(floodst)

Question 6:

plot_rsum = sum(floodst)

raster_stats = cellStats(floodst, stat=sum)*res(floodst)^2/1e6 #showing area covered in km&2

Step 1

kable(raster_stats, col.names = c("Area Covered"), caption = "Area Covered by each classification") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Area Covered by each classification
Area Covered
layer.1 5.9994
layer.2 6.4908
layer.3 10.7451
layer.4 7.6221
layer.5 13.6809
layer.6 7.9371

Step 2

plot(plot_rsum, col = blues9)

Step 3

plot_copy = (plot_rsum[plot_rsum == 0] = NA)
plot(plot_rsum, col = blues9)

Why is our raster STACK, comprised of 6 RASTERS, have values ranging from 6 to 0?

  • Perhaps they’re an average of the values at that pixel.

Extra Credit:

Used mapview from previous question to identify area of interest

aoi = st_point(c(-91.78946, 42.06303)) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf() %>% 
  st_transform(st_crs(floodst))

extract(floodst, aoi)
##      layer.1 layer.2 layer.3 layer.4 layer.5 layer.6
## [1,]       0       0       1       0       1       1

Here I identified which layers correctly labeled the location in the image as a flood.