Geography 176A

Lab 05: Rasters and Remote Sensing

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

# Data
library(USAboundaries)
library(rnaturalearth)

# Visualization
library(gghighlight)
library(ggrepel)
library(knitr)
library(rmapshaper)
library(mapview)
library(leaflet)
library(elevatr)
library(raster)
library(getlandsat)

knitr::opts_chunk$set(cache = F)

Question 1

setwd("C:/Users/slaye/Documents/github/geog-176A-labs/")
bb=read_csv("data/uscities.csv") %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>% 
  filter(city == "Palo", state_name == "Iowa") %>% 
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()

Question 2

#2.2
setwd("C:/Users/slaye/Documents/github/geog-176A-labs/")
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)

#2.3
st = sapply(files, lsat_image)

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

#2.4
cropper = bb %>% st_transform(crs(s))

r = crop(s, cropper)

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      : band1, band2, band3, band4, band5, band6 
## 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) crs : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs resolution : 30, 30 (x, y)

Question 3

#1.1
par(mfrow = c(2,2))
RGB = plotRGB(r, r=4, g=3, b=2)
NIRRG = plotRGB(r, r=5, g=4, b=3)
NIRSWIR1G = plotRGB(r, r=5, g=6, b=4)
TIRSCoastalCirrus = plotRGB(r, r=10, g=1, b=9)

#3.2
par(mfrow = c(2,2))
RGB2 = plotRGB(r, r=4, g=3, b=2, stretch = "lin")
NIRRG2 = plotRGB(r, r=5, g=4, b=3, stretch = "lin")
NIRSWIR1G2 = plotRGB(r, r=5, g=6, b=4, stretch = "hist")
TIRSCoastalCirrus2 = plotRGB(r, r=10, g=1, b=9, stretch = "hist")

The stretch improves the contrast in th emiage so you can see the color differences better.

Question 4

#4.1
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)

palette = colorRampPalette(c("blue", "white", "red"))

(water = stack(ndvi, ndwi, mndwi, wri, swi) %>% 
  setNames(c("vegetation", "water", "modwater", "waterratio", "simplewater")))

## class      : RasterStack 
## dimensions : 340, 346, 117640, 5  (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 
## names      :  vegetation,       water,    modwater,  waterratio, simplewater 
## min values : -0.12413047, -0.51084642, -0.33825365,  0.39078526,  0.01635175 
## max values :   0.5779096,   0.1793064,   0.2033177,   1.3475451,         Inf

plot(water, col = palette(256))

#4.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)}

f1 = calc(ndvi, thresholding1)
f2 = calc(ndwi, thresholding2)
f3 = calc(mndwi, thresholding3)
f4 = calc(wri, thresholding4)
f5 = calc(swi, thresholding5)

(flood = stack(f1, f2, f3, f4, f5) %>% 
  setNames(c("vegetation", "water", "modwater", "waterratio", "simplewater")))

## class      : RasterStack 
## dimensions : 340, 346, 117640, 5  (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 
## names      : vegetation, water, modwater, waterratio, simplewater 
## min values :          1,     1,        1,          1,           1 
## max values :          1,     1,        1,          1,           1

plot(flood, col = "blue")

Question 5

#5.1
set.seed

## function (seed, kind = NULL, normal.kind = NULL, sample.kind = NULL) 
## {
##     kinds <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper", 
##         "Mersenne-Twister", "Knuth-TAOCP", "user-supplied", "Knuth-TAOCP-2002", 
##         "L'Ecuyer-CMRG", "default")
##     n.kinds <- c("Buggy Kinderman-Ramage", "Ahrens-Dieter", "Box-Muller", 
##         "user-supplied", "Inversion", "Kinderman-Ramage", "default")
##     s.kinds <- c("Rounding", "Rejection", "default")
##     if (length(kind)) {
##         if (!is.character(kind) || length(kind) > 1L) 
##             stop("'kind' must be a character string of length 1 (RNG to be used).")
##         if (is.na(i.knd <- pmatch(kind, kinds) - 1L)) 
##             stop(gettextf("'%s' is not a valid abbreviation of an RNG", 
##                 kind), domain = NA)
##         if (i.knd == length(kinds) - 1L) 
##             i.knd <- -1L
##     }
##     else i.knd <- NULL
##     if (!is.null(normal.kind)) {
##         if (!is.character(normal.kind) || length(normal.kind) != 
##             1L) 
##             stop("'normal.kind' must be a character string of length 1")
##         normal.kind <- pmatch(normal.kind, n.kinds) - 1L
##         if (is.na(normal.kind)) 
##             stop(gettextf("'%s' is not a valid choice", normal.kind), 
##                 domain = NA)
##         if (normal.kind == 0L) 
##             stop("buggy version of Kinderman-Ramage generator is not allowed", 
##                 domain = NA)
##         if (normal.kind == length(n.kinds) - 1L) 
##             normal.kind <- -1L
##     }
##     if (!is.null(sample.kind)) {
##         if (!is.character(sample.kind) || length(sample.kind) != 
##             1L) 
##             stop("'sample.kind' must be a character string of length 1")
##         sample.kind <- pmatch(sample.kind, s.kinds) - 1L
##         if (is.na(sample.kind)) 
##             stop(gettextf("'%s' is not a valid choice", sample.kind), 
##                 domain = NA)
##         if (sample.kind == 0L) 
##             warning("non-uniform 'Rounding' sampler used", domain = NA)
##         if (sample.kind == length(s.kinds) - 1L) 
##             sample.kind <- -1L
##     }
##     .Internal(set.seed(seed, i.knd, normal.kind, sample.kind))
## }
## <bytecode: 0x000000003e0d8b88>
## <environment: namespace:base>

#5.2
rvalues = getValues(r) %>% 
  na.omit
dim(rvalues)

## [1] 117640      6

#these values indicate that the values were extracted as points from 2 separate groups
par(mfrow = c(2,2))
kave = stats::kmeans(rvalues, 4)
kmeans_raster = flood$water
values(kmeans_raster) = kave$cluster
plot(kmeans_raster)

kave = stats::kmeans(rvalues, 8)
kmeans_raster = flood$water
values(kmeans_raster) = kave$cluster
plot(kmeans_raster)

kave = stats::kmeans(rvalues, 12)
kmeans_raster = flood$water
values(kmeans_raster) = kave$cluster
plot(kmeans_raster)

kave = stats::kmeans(rvalues, 16)
kmeans_raster = flood$water
values(kmeans_raster) = kave$cluster
plot(kmeans_raster)

#5.3
kave = stats::kmeans(rvalues, 14)
kmeans_raster = flood$water
values(kmeans_raster) = kave$cluster
plot(kmeans_raster)

f1[is.na(f2[])] <- 0 
f2[is.na(f2[])] <- 0 
f3[is.na(f2[])] <- 0 
f4[is.na(f2[])] <- 0 
f5[is.na(f2[])] <- 0 

floodtable = table(values(f2), values(kmeans_raster))
floodtable

##    
##         1     2     3     4     5     6     7     8     9    10    11    12
##   0  8001  7640  1034  4001 22424  7610   660  3350  4941  7118  5560  8427
##   1     0     0  7156     0     0     0     1     0     0     0    52     0
##    
##        13    14
##   0 20466  9197
##   1     2     0

which.max(floodtable)

## [1] 9

thresholding6 = function(x){ifelse(x > 5, 1, NA)}
f6 = calc(kmeans_raster, thresholding6)

flood2 = stack(f1, f2, f3, f4, f5, f6) %>% 
  setNames(c("vegetation", "water", "modwater", "waterratio", "simplewater", "kmeans"))
plot(flood2, col = "blue")

Question 6

stats = cellStats(flood2, sum)
stats = (stats*900)/1000
knitr::kable(stats, 
             caption = "Area in km^2 Affected by Flooding per Raster",
            
             format.args = list(big.mark = ","))

Area in km^2 Affected by Flooding per Raster

	x
vegetation	5,985.9
water	6,489.9
modwater	10,742.4
waterratio	7,621.2
simplewater	13,680.9
kmeans	60,645.6

uncertain = calc(flood2, fun = sum)

plot(uncertain)

mapview(uncertain)

Some of the cell values are not an even number because the data is averaged over the raster box, meaning the boxes were partially overlapped somehow which leads to a number that is not an integer.

Extra Credit

bb2 = st_transform(bb, 4326)

osm = osmdata::opq(bb2) %>% 
  osmdata::add_osm_feature("building") %>% 
  osmdata::osmdata_sf()
leaflet() %>% 
  setView(lng=-91.78943, lat=42.06307, zoom = 13) %>% 
  addProviderTiles(providers$CartoDB) 

local <- st_point(c(-91.78943, 42.06307))

sf_local = st_sfc(local, crs = st_crs(flood2)) %>% 
  st_cast("POINT")

sf_local

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -91.78943 ymin: 42.06307 xmax: -91.78943 ymax: 42.06307
## CRS:            +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs

floodvalue = raster::extract(flood2, sf_local, na.rm = TRUE)

## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'extract' for signature '"RasterStack", "sfc_POINT"'

flood2

## class      : RasterStack 
## 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 
## names      : vegetation, water, modwater, waterratio, simplewater, kmeans 
## min values :          0,     0,        1,          1,           1,      1 
## max values :          1,     1,        1,          1,           1,      1