library(sf)        # vector manipulation
library(raster)    # raster manipulation
library(whitebox)  # terrain analysis
library(osmdata)   # OSM API
library(elevatr)   # Elevation  Web Tiles
library(tidyverse)
library(units)
library(fasterize)

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

elev = elevatr::get_elev_raster(bound, z = 13) %>% 
  crop(bound)  

writeRaster(elev, "../data/area-elev.tif", overwrite = TRUE)

el = raster("../data/area-elev.tif")

plot(el)

# Buildings and river-network data

# building

bb = st_bbox(bound) %>% 
  st_as_sfc() %>% 
  st_as_sf() %>% 
  st_transform(4326) 

build = osmdata::opq(bb) %>% 
  add_osm_feature(key = 'building') %>% 
  osmdata_sf()

bd = build$osm_polygons %>% 
   dplyr::select(osm_id, name, building) 

cent = st_centroid(bd) %>% 
  st_transform(crs(bound))

cent_bound = st_intersection(cent, bound)

# stream

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

poly = stream$osm_points %>% 
  dplyr::select(osm_id) %>% 
  st_combine() %>% 
  st_cast('POLYGON')

str = stream$osm_lines %>% 
   dplyr::select(osm_id, name, waterway)

str_bound = st_intersection(str, bb)

plot(str_bound$geometry)

# Terrain Analysis

# Hillshade

wbt_hillshade("../data/area-elev.tif", "../data/area-hillshade.tif")
## [1] "hillshade - Elapsed Time (excluding I/O): 0.38s"
r = raster("../data/area-hillshade.tif")

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

plot(bound)
plot(str_bound$geometry, add = TRUE, col = 'darkred')

# Height Above Nearest Drainage

# Creating the river raster