pretty_rivers/02_process_data.R

44 lines
1.3 KiB
R
Raw Normal View History

# load libraries
2023-10-16 16:46:49 -05:00
library(nhdplusTools)
library(sf)
2023-10-16 15:24:53 -05:00
# 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(
2023-10-16 16:59:19 -05:00
x = list(
cbind(
extent$longitude[c(1,2,2,1,1)],
extent$latitude[c(1,1,2,2,1)])
)
2023-10-16 15:24:53 -05:00
)
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),
2023-10-16 16:59:19 -05:00
longitude_min = min(extent$longitude),
latitude_max = max(extent$latitude),
latitude_min = min(extent$latitude))
2023-10-16 15:24:53 -05:00
# load and crop hydrologic data
data <- list(NULL)
sf_use_s2(FALSE)
2023-10-16 16:46:49 -05:00
layers <- c("NHDArea",
2023-10-16 16:59:19 -05:00
"NHDFlowline",
"NHDWaterbody",
"NHDPlusLandSea")
for (layer in layers){
2023-10-16 16:59:19 -05:00
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) {
2023-10-16 16:59:19 -05:00
political[[boundary]] <- st_read(paste0(data_dir, "/political_boundaries/", boundary))
}