#Library Codes
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(osmdata)
library(kableExtra)

#Question 1:

bb = read_csv("~/github/geog176A-daily-exercises/data/uscities.csv") %>% #reading in csv file 
  filter(city == "Palo")  %>% 
  #Filtering city to Palo, Iowa 
  st_as_sf(coords = c("lng", "lat"), crs = 4326)%>% 
  #Structuring to spatial features by defining coordinate fields and CRS
  st_transform(5070)  %>% 
  st_buffer(5000)  %>% #Generating 5 kilometer buffer around point
  st_bbox()  %>% #finding the bounding box of the buffered region
  st_as_sfc() %>% 
  st_as_sf()

# mapview(bb)
bbwgs = bb %>% 
  st_transform(4326)

bb = st_bbox(bbwgs)

#Question 2

#2.1 loading all available scenes and assigning it to an object
scenes = lsat_scenes()
#Filerting loaded data to include only those that are suitable for our AOI on the data 2016-09-26
down = scenes %>% 
  filter(min_lat <= bb$ymin, max_lat >= bb$ymax, min_lon <= bb$xmin, max_lon >= bb$xmax, as.Date(acquisitionDate) == as.Date("2016-09-26"))
#saving data to a csv file
write.csv(down, file = "/Users/Allan/github/geog-176A-labs/data/palo-flood.csv", row.names = FALSE)
#2.2 Downloading cached data on our computer
files = lsat_scene_files("https://s3-us-west-2.amazonaws.com/landsat-pds/L8/025/031/LC80250312016270LGN00/index.html") %>% 
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file))  %>% #passing downloaded url
  arrange(file) %>% 
  pull(file)
#2.3 apply last image function over our URL path to return a vector of local file paths
st = sapply(files, lsat_image)
#creating raster stack
s = stack(st) %>% 
  setNames(c(paste0("band", 1:6)))

The dimension of our stacked image is 7811 roqs x 7681 columns with 59996291 cells and 6 layers. The crs is in +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs. The resolution is 30, 30 (x, y).

#2.4 cropping raster stack to analyze our image for the regions surroud AOI
cropper = bbwgs %>% 
  st_transform(crs(s))
r = crop(s, cropper)
dim(r)
## [1] 340 346   6

The dimensions of our cropped image stack, r, are 340 row x 346 columns with 117640 cells and 6 layers. The crs is in “+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs”. The resolution is 30, 30 (x, y).

#Question 3

#3.1 Working with Landsat bands to make combination of color plots for our AOI

#3.2 Replicating images by applying stretch to improve clarity and/or contrast
par(mfrow = c(2,1))
#R-G-B with Lin/Hist Stretches
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
plotRGB(r, r = 4, g = 3, b = 2, stretch = "hist")

#NIR-R-G with Lin/Hist Stretches
plotRGB(r, r = 5, g = 5, b = 3, stretch = "lin")
plotRGB(r, r = 5, g = 5, b = 3, stretch = "hist")

#NIR-SWIR1-R with Lin/Hist Stretches
plotRGB(r, r = 5, g = 6, b = 4, stretch = "lin")
plotRGB(r, r = 5, g = 6, b = 4, stretch = "hist")

#NIR-G-B with Lin/Hist Stretches
plotRGB(r, r = 5, g = 3, b = 2, stretch = "lin")
plotRGB(r, r = 5, g = 3, b = 2, stretch = "hist")

dev.off()
## null device 
##           1

Application of the color stretch is used to change the clarity/contrast of our images. These 2 images are then compared to one another to highlight if map features are more noticable on one image than the other, and what these differences highlight.

#Question 4

#4.1
#Making our palette for raster stacking. 
palette = colorRampPalette(c("blue", "white", "red"))

#creating rasters using our formulas
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))

stack_indices = raster::stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(stack_indices)

### The 5 images are all the same but are being plotted using different color bands. They are meant to highlight the flood region that we are trying to analyze. They highlight the area of interest using different pigments. They deviate in how well the colors being shown can bring out the flood region better. in the ndvi image, its hard to see the flood region at all, but in our other 4 images, the flood region is outlined pretty evidently. The 5th image uses shades of white to outline the flood region.

#4.2 Creating Raster Thresholding
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)}
#Using calc functions to apply a custom formula using the thresholdings above
flood1 = calc(ndvi, thresholding1)
plot(flood1, col = "blue")

flood2 = calc(ndwi, thresholding2)
plot(flood2, col = "blue")

flood3 = calc(mndwi, thresholding3)
plot(flood3, col = "blue")

flood4 = calc(wri, thresholding4)
plot(flood4, col = "blue")

flood5 = calc(swi, thresholding5)
plot(flood5, col = "blue")

#stacking the binary files into a new stack.
floods = raster::stack(flood1, flood2, flood3, flood4, flood5) %>% 
  setNames(c("ndvi", "ndwi", "mndwi", "wri", "swi"))

#plotting 
plot(floods, col = blues9)

#Question 5

#5.1
#Setting seed.
set.seed(09062020)

#5.2 Extracting values from our 6-banded raster stack
r_values = getValues(r) %>% 
  na.omit()

#clustering extracted data to a given numer of centers, k. 
k_12 = kmeans(r_values, 12)
k_6 = kmeans(r_values, 6)
k_4 = kmeans(r_values, 4)

#Creating new raster object by copying one of the original bands. 
kmeans_raster = stack_indices$NDVI
kmeans_raster6 = stack_indices$NDVI
kmeans_raster4 = stack_indices$NDVI

#setting the values of the copied raster to the cluster vector from the output kmeans object. 
values(kmeans_raster) = k_12$cluster
values(kmeans_raster6) = k_6$cluster
values(kmeans_raster4) = k_4$cluster

#Plotting.
plot(kmeans_raster)

plot(kmeans_raster6)

plot(kmeans_raster4)

### The dimensions of the extracted values tells us that the data was extracted by pulling the numnerical integer value of each cell and turned them all into clusters by band number.

#5.3 
#Generating a table crossing the values of one of our binary flood rasters
floods_values = values(flood1)
k_12_values = values(kmeans_raster)
kmeans_table = table(floods_values, k_12_values)

idx = which.max(kmeans_table)

kmeans_raster[kmeans_raster != idx] = 0

kmeans_stack = stack(flood1, flood2, flood3, flood4, flood5, kmeans_raster)

plot(kmeans_stack)

#Question 6

#Making a sum function of our last stack from Q5
plot_sum = calc(kmeans_stack, fun = sum)

floods_stats = cellStats(kmeans_stack, stat = sum)*res(kmeans_stack)^2/1e6

kable(floods_stats, col.names = "Area (m^2)")
Area (m^2)
layer.1 5.9994
layer.2 6.4908
layer.3 10.7451
layer.4 7.6221
layer.5 13.6809
NDVI 49.1670

Some of the cell values are not an even number because the NA values have been removed from all of these cells, and thus, some bands might have lost an odd number of cells as a result.

plot(plot_sum, col = blues9)

plot_sum[plot_sum == 0] = NA
mapview(plot_sum)

#Extra Credit

new_aoi1 = st_point(c(-91.78946, 42.06303)) %>% 
  st_sfc(crs = 4326) %>% 
  st_as_sf() %>% 
  st_transform(st_crs(kmeans_stack))

new_aoi2 = st_point(c(-91.78482, 42.06245)) %>% 
  st_sfc(crs = 4326) %>% 
  st_as_sf() %>% 
  st_transform(st_crs(kmeans_stack))

f_values = raster::extract(kmeans_stack, new_aoi1) 
f_values_plot = new_aoi1 %>% 
  mutate(f_values = f_values)