Collecting Data

basin = read_sf("https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-11119750/basin/")
elev = elevatr::get_elev_raster(basin, z = 13) %>%
  crop(basin) %>% 
  mask(basin)
 
writeRaster(elev, filename =  "../data/basin_elev.tif", overwrite = TRUE)
new_units = elev * 3.281

osm = osmdata::opq(basin)%>% 
  add_osm_feature(key= 'building') %>% 
  osmdata_sf()
 
polygon = osm$osm_polygons %>% 
   st_centroid()

clip_build = st_intersection(polygon, basin)

railway = clip_build %>% 
  dplyr::filter(amenity == "railway")

streams = osmdata::opq(basin)%>% 
  add_osm_feature(key= 'waterway', value= "stream") %>% 
  osmdata_sf()

line = streams$osm_lines 

clip_streams = st_intersection(line, basin)

Terrain Analysis

wbt_hillshade("../data/basin_elev.tif", "../data/basin_hillshade.tif")
## [1] "hillshade - Elapsed Time (excluding I/O): 0.32s"
hillshade = raster( "../data/basin_hillshade.tif")

plot(hillshade, col= gray.colors(256, alpha = .5), legend = FALSE)

plot(basin, add= TRUE, lwd = 2)

plot(clip_streams, add= TRUE, lwd = 2, col = "blue")

stream_buff = clip_streams %>% 
  st_transform(5070) %>% 
  st_buffer(10) %>% 
  st_transform(4326)
  
stream_fast = fasterize::fasterize(stream_buff, new_units) 

writeRaster(stream_fast,"../data/stream_raster.tif", overwrite = TRUE)
wbt_breach_depressions("../data/basin_elev.tif", "../data/basin_breach_depressions.tif")
## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.161s"
wbt_elevation_above_stream("../data/basin_breach_depressions.tif", "../data/stream_raster.tif",  "../data/hand.tif")
## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.53s"

Correcting to local reference datum

hand_raster = raster("../data/hand.tif")

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

new_hand =  hand_raster + 3.69

new_hand[stream_fast == 1] = 0

writeRaster(new_hand,"../data/hand_offset.tif", overwrite = TRUE)

2017 Impact Assessment

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

hand_offset[hand_offset > 10.02] = NA


plot(hillshade, col= gray.colors(256, alpha = .5), legend = FALSE, axes = FALSE, box = FALSE )

plot(hand_offset, add = TRUE, col = rev(blues9))

plot(railway, add = TRUE, col = "green", cex = 1, pch = 16)

building_cen =  ifelse(!is.na(raster::extract(hand_offset, polygon)), "red", "black")


plot(hillshade, axes = FALSE, box = FALSE, col= gray.colors(256, alpha = .5), legend = FALSE, main = paste(sum(building_cen == "red"),  "Impacted Buildings"))

plot(hand_offset, add = TRUE, col = rev(blues9))

plot(clip_build$geometry, add = TRUE, col = building_cen, pch = 16, cex = .08)

plot(railway, add= TRUE, col= "green", cex = 1, pch = 16)