# load libraries library(nhdplusTools) library(sf) # reload extent data so you don't need to run the download data script again data_dir <- "data" nhdplusTools_data_dir(dir = paste0(data_dir, "/hydrologic")) extent <- read.csv(file = "extent.csv") crs <- 4269 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) extent_bbox <- st_bbox(st_transform(x = extent_poly, crs)) extent <- list(longitude_max = max(extent$longitude), longitude_min = min(extent$longitude), latitude_max = max(extent$latitude), latitude_min = min(extent$latitude)) # 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)