44 lines
1.2 KiB
R
44 lines
1.2 KiB
R
# 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) |