# load libraries
library(nhdplusTools)
library(sf)

# set hydrologic data dir
data_dir <- "data"
ifelse(
  !dir.exists(file.path(getwd(), paste0(data_dir, "/hydrologic"))),
  dir.create(file.path(getwd(), paste0(data_dir, "/hydrologic"), recursive = TRUE)), 
  FALSE)

nhdplusTools_data_dir(dir = paste0(data_dir, "/hydrologic"))

# load extent of map
extent <- read.csv(file = "extent.csv")
crs <- 4269

# make extent a polygon
extent_poly <- st_polygon(
  x = list(
    cbind(
      extent$longitude[c(1,2,2,1,1)], 
      extent$latitude[c(1,1,2,2,1)])
  )
)

extent_poly <- st_sfc(extent_poly, crs=4326)
# make extent a bounding box
extent_bbox <- st_bbox(st_transform(x = extent_poly, crs))

# make list of extent coordinates
extent <- list(longitude_max = max(extent$longitude),
               longitude_min = min(extent$longitude),
               latitude_max = max(extent$latitude),
               latitude_min = min(extent$latitude))