# load libraries
library(nhdplusTools)
library(sf)

# get watershed areas
extent_huc <- get_huc(AOI = extent_poly, buffer = 0, type = "huc04")

# download data
download_nhdplushr(nhdplusTools_data_dir(), extent_huc$huc4, download_files = TRUE)

# get vaa
get_vaa()

# combine data into geopackage
layers <- c("NHDArea",
            "NHDFlowline",
            "NHDWaterbody",
            "NHDPlusLandSea")

data <- get_nhdplushr(hr_dir = nhdplusTools_data_dir(), 
                      out_gpkg = paste0(nhdplusTools_data_dir(),"/data.gpkg"),
                      layers = layers,
                      overwrite = TRUE)