library(sf)        # vector manipulation
library(raster)    # raster manipulation
library(fasterize) # "faster" raster
library(whitebox)  # terrain analysis

# Data libraries
library(osmdata)   # OSM API
library(elevatr)   # Elevation  Web Tiles

library(tidyverse)
library(units)
library(knitr)

#Question1

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

#Question 2

wbt_hillshade("../data/basin-elev.tif", "../data/basin-hillshade.tif")
## [1] "hillshade - Elapsed Time (excluding I/O): 0.28s"
hill_r = raster("../data/basin-hillshade.tif")
plot(hill_r, axes = FALSE, box = FALSE, col = gray.colors(256, alpha = 0.5), main = "Hillshade")
plot(stream, add = TRUE, col = "blue")

#Question2

## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.290s"
## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.35s"
plot(hill_r, axes = FALSE, box = FALSE, col = gray.colors(256, alpha = 0.5), legend = FALSE)
plot(HAND_offset, add = TRUE, col = rev(blues9))
plot(railway, add = TRUE, col = "green", cex = 1, pch = 16)