Lab 06

library(sf)
library(raster)
library(fasterize)
library(whitebox)
library(osmdata)
library(readr)
library(elevatr)
library(mapview)
library(dplyr)
library(gifski)

Collecting Data

Basin boundary

basin_boundary = read_sf("https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-11119750/basin")

Elevation Data

basin_elev = elevatr::get_elev_raster(basin_boundary, z = 13) %>%
  crop(basin_boundary) %>% 
  mask(basin_boundary) 

basin_feet = basin_elev *3.281

writeRaster(basin_feet, filename = "C:/Users/hopew/Desktop/github176/geog-176A-labs/data/basin_elevation_ft.tif", overwrite = T)

Buildings and river-network data

building = osmdata::add_osm_feature(opq(basin_boundary), "building") %>% 
  osmdata_sf()

stream = osmdata::add_osm_feature(opq(basin_boundary), "waterway", "stream") %>% osmdata_sf()

building = st_centroid(building$osm_polygons) %>% 
  st_intersection(basin_boundary)

railway = dplyr::filter(building, amenity =="railway")

river = st_intersection(stream$osm_lines, basin_boundary)

Terrain Analysis

Hillshade

wbt_hillshade("C:/Users/hopew/Desktop/github176/geog-176A-labs/data/basin_elevation_ft.tif", "C:/Users/hopew/Desktop/github176/geog-176A-labs/data/basin_hillshade.tif")

## [1] "hillshade - Elapsed Time (excluding I/O): 0.75s"

basin_hillshade = raster("C:/Users/hopew/Desktop/github176/geog-176A-labs/data/basin_hillshade.tif")

plot(basin_hillshade, col = gray.colors(256, alpha = .5), legend = F)
plot(river$geometry, col = "blue", add = T)
plot(basin_boundary$geometry, add = T)

Height Above Nearest Drainage

Creating the river raster

river_raster = st_transform(river, 5070) %>% 
  st_buffer(10) %>% 
  st_transform(crs(basin_feet)) %>% 
  fasterize::fasterize(basin_feet) %>% 
  writeRaster("../data/river_raster.tif", overwrite = T)

Creating the hydrologically corrected surface

wbt_breach_depressions("../data/basin_elevation_ft.tif",
                       "../data/breach_depressions.tif")

## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.309s"

Creating the HAND raster

wbt_elevation_above_stream("../data/breach_depressions.tif", "../data/river_raster.tif", "../data/elev_above_stream.tif")

## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.89s"

Correcting to local reference datum

hand_raster = raster("../data/elev_above_stream.tif") +3.69

river_raster = raster("../data/river_raster.tif")

hand_raster[river_raster ==1] = 0

writeRaster(hand_raster, "../data/corrected.tif", overwrite = T)

2017 Impact Assessment

Map the flood and estimate the impacts

flood17 = raster("../data/corrected.tif")
flood17_map = flood17

flood17_map[flood17_map >= 10.02] = NA

cols = ifelse(!is.na(raster::extract(flood17_map, building)), "red", "black")


plot(basin_hillshade, col = gray.colors(256, alpha = .5), legend = F, main = paste0(sum(cols == "red"), " impacted buildings, 10.02 foot stage."))

plot(flood17_map, col = rev(blues9), legend = F, add = T)
plot(building$geometry, add = T, col = cols, pch = 16, cex = .08)



plot(railway$geometry, col = "green", cex = 1, pch = 16, add = T)

plot(basin_boundary$geometry, add = T, border = "black")

Does the map look accurate?

The map looks accurate.

Extra Credit

Flood Inundation Map library

sb = AOI::aoi_get("Santa Barbara")
flood_gif = crop(flood17, sb)
basin_clip = st_intersection(basin_boundary, sb)
hill_crop = crop(basin_hillshade, sb)

Making the GIF

gifski::save_gif({
  for(i in 0:20){
    tmp = flood_gif
    tmp[tmp>= i] = NA
    cols = ifelse(!is.na(raster::extract(tmp, building)), "red", "black")
    plot(hill_crop, col = gray.colors(256, alpha = .5), legend = F,
         main = paste0(sum(cols == "red"), " impacted buildings, ", i, " foot stage"))
    plot(tmp, col = rev(blues9), legend = F, add = T)
    plot(building$geometry, add = T, col = cols, pch = 16, cex = .08)
    plot(railway$geometry, add = T, col = "green", cex = 1, pch = 16)
    plot(basin_clip$geometry, add = T, border = "black")
  }
}, gif_file = "../data/mission_creek.gif",
width = 600, height = 600,
delay = .7, loop = T
)

## 
Frame 1 (4%)
Frame 2 (9%)
Frame 3 (14%)
Frame 4 (19%)
Frame 5 (23%)
Frame 6 (28%)
Frame 7 (33%)
Frame 8 (38%)
Frame 9 (42%)
Frame 10 (47%)
Frame 11 (52%)
Frame 12 (57%)
Frame 13 (61%)
Frame 14 (66%)
Frame 15 (71%)
Frame 16 (76%)
Frame 17 (80%)
Frame 18 (85%)
Frame 19 (90%)
Frame 20 (95%)
Frame 21 (100%)
## Finalizing encoding... done!

## [1] "C:\\Users\\hopew\\Desktop\\github176\\geog-176A-labs\\data\\mission_creek.gif"