Libraries

library(sf)
library(raster)
library(fasterize)
library(whitebox)
library(osmdata)
library(elevatr)

Data Gathering

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

Elevation Data

elev  = get_elev_raster(basin, z = 13) %>%
  crop(basin) %>%
  mask(basin)

elev_ft = elev * 3.281

writeRaster(elev_ft, filename = "../data/goleta-creek-elev.tif", overwrite = TRUE)

Buildings and river-network data

q = opq(basin)
building =  add_osm_feature(q, "building") %>% osmdata_sf()
water = add_osm_feature(q, "waterway", "stream")  %>% osmdata_sf()

build = st_centroid(building$osm_polygons) %>% st_intersection(basin)

railroad = dplyr::filter(build, amenity == "railway")

riv =  st_intersection(water$osm_lines, basin)

Terrain Analysis

Hillshade

wbt_hillshade("../data/goleta-creek-elev.tif", 
              "../data/goleta-creek-hill.tif",
              zfactor = 0.00000316)
## [1] "hillshade - Elapsed Time (excluding I/O): 0.11s"
hillshade = raster("../data/goleta-creek-hill.tif")

par(mar=c(0,0,0,0))
plot(hillshade, col = gray.colors(256, alpha = .5), legend = FALSE, box = FALSE, axes = FALSE)
plot(basin$geometry, add = TRUE)
plot(riv$geometry, col = "navy", add = TRUE)

Creating the river raster

utf = st_transform(riv, 5070) %>%
  st_buffer(10) %>%
  st_transform(crs(elev_ft)) %>% 
  fasterize(elev_ft)

writeRaster(utf, filename = "../data/river.tif", overwrite = TRUE)

Creating the hydrologically corrected surface

wbt_breach_depressions("../data/goleta-creek-elev.tif",  "../data/goleta-creek-b.tif")
## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.116s"

Creating the HAND raster

wbt_elevation_above_stream("../data/goleta-creek-b.tif", "../data/river.tif", '../data/hand.tif')
## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.30s"

Correcting to local reference datum

h = raster('../data/hand.tif')
r = raster('../data/river.tif')
h = h + 3.69
h[r == 1] = 0
writeRaster(h, "../data/hand-corrected.tif", overwrite = TRUE)

2017 Impact Assessment

Mapping the flood

Estimate the impacts

h = raster("../data/hand-corrected.tif")
tmp = h
tmp[tmp >= 10.02] = NA
cols = ifelse(!is.na(raster::extract(tmp, build)), "red", "black")

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

    plot(tmp, col = rev(blues9), legend = FALSE, add = TRUE)
    
    plot(build$geometry, add = TRUE, pch = 16, cex = .06, col = cols)
    
    plot(railroad$geometry, col = "green", pch = 16, add = TRUE)
    
    plot(basin$geometry, add = TRUE, border = "black")

Flood Inudation Map library

(sb = AOI::aoi_get("Santa Barbara"))
## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -119.8598 ymin: 34.33603 xmax: -119.6399 ymax: 34.46392
## geographic CRS: WGS 84
##                         geometry       request
## 1 POLYGON ((-119.8598 34.3360... Santa Barbara
h = raster("../data/hand-corrected.tif")
bas2 = st_intersection(basin, sb)
hill2 = crop(hillshade, sb)
h2 = crop(h, sb)

gifski::save_gif({
  for(t in 0:20) {
    tmp = h2
    tmp[tmp > t] = NA
    cols = ifelse(!is.na(raster::extract(tmp, build)), "red", "black")
    plot(hill2, col = gray.colors(256, alpha = .5),  legend = FALSE, 
         main = paste0(sum(cols == "red"), " impacted buildings, ", t, " foot stage"))
    plot(tmp, col = rev(blues9), legend = FALSE, add = TRUE)
    plot(build$geometry, add = TRUE, pch = 16, cex = .25, col = cols)
    plot(railroad$geometry, col = "green", pch = 16, add = TRUE)
    plot(bas2$geometry, add = TRUE, border = "black")
  }
}, gif_file = "../data/cool.gif", width = 600, height = 600, delay = .7, loop = TRUE)
## 
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\\Ryan Erickson\\Documents\\GitHub\\geog176A_labs\\data\\cool.gif"

../data/cool.gif