# load libraries
library(nhdplusTools)
library(sf)
# load and crop hydrologic data
data <- list(NULL)
sf_use_s2(FALSE)
layers <- c("NHDArea",
            "NHDFlowline",
            "NHDWaterbody",
            "NHDPlusLandSea")
for (layer in layers){
  data[[layer]] <- st_crop(st_read(paste0(nhdplusTools_data_dir(),"/data.gpkg"), layer = layer), y = extent_bbox)
}

# load political boundaries
political <- st_crop(
  st_transform(
    st_read(
      paste0(data_dir, "/political_boundaries/bound_p/boundaries_p_2021_v3.shp")),
    crs=crs),
  y = extent_bbox)