# load libraries 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) 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 <- list(NULL) political_boundaries <- list.files(path = paste0(data_dir, "/political_boundaries")) for (boundary in political_boundaries) { political[[boundary]] <- st_read(paste0(data_dir, "/political_boundaries/", boundary)) }