initial 02_process_data.R and 03_make_figures.R

This commit is contained in:
Ben Varick 2023-10-16 12:17:24 -05:00
parent b0a038504f
commit 7eddee51b4
3 changed files with 72 additions and 3 deletions

View File

@ -1,11 +1,12 @@
# load libraries
library(tidyverse)
library(nhdplusTools)
library(sf)
nhdplusTools_data_dir(dir = "data")
# load extent of map
extent <- read.csv(file = "extent.csv")
crs <- 26916
crs <- 4269
extent_poly <- st_polygon(
x = list(
@ -16,7 +17,7 @@ extent_poly <- st_polygon(
)
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),
@ -28,3 +29,8 @@ extent_huc <- get_huc(AOI = extent_poly, buffer = 0, type = "huc04")
# download data
download_nhdplushr("data", extent_huc$huc4, download_files = TRUE)
# combine data into geopackage
data <- get_nhdplushr(hr_dir = "data",
out_gpkg = "data/data.gpkg",
layers = NULL)

View File

@ -0,0 +1,16 @@
# load libraries
library(sf)
# load data
layers <- c("NHDArea",
"NHDFlowline",
"NHDWaterbody",
"NHDPlusLandSea")
data <- list(NULL)
sf_use_s2(FALSE)
for (layer in layers){
data[[layer]] <- st_crop(st_read("data/data.gpkg", layer = layer), y = extent_bbox)
}

View File

@ -0,0 +1,47 @@
# load libraries
library(sf)
colors <- list(darkblue = "#062e57",
lightblue = "#b1dcf3")
# save figure
tiff(filename = "figures/map.tiff",
width = 6600,
height = 10200,
res = 600,
compression = "lzw")
# plot map
plot(sf::st_geometry(extent_poly),
col = colors$darkblue,
extent = extent_bbox,
xlim = c(extent$longitude_min, extent$longitude_max),
ylim = c(extent$latitude_min, extent$latitude_max),
border = "black")
plot(sf::st_geometry(data$NHDWaterbody),
col = colors$lightblue,
border = NA,
extent = extent_bbox,
xlim = c(extent$longitude_min, extent$longitude_max),
ylim = c(extent$latitude_min, extent$latitude_max),
add = TRUE)
plot(sf::st_geometry(data$NHDArea),
col = colors$lightblue,
border = NA,
extent = extent_bbox,
xlim = c(extent$longitude_min, extent$longitude_max),
ylim = c(extent$latitude_min, extent$latitude_max),
add = TRUE)
plot(sf::st_geometry(data$NHDFlowline),
col = colors$lightblue,
lwd = data$NHDFlowline$TotDASqKM^0.3204*0.0446,
border = NA,
extent = extent_bbox,
xlim = c(extent$longitude_min, extent$longitude_max),
ylim = c(extent$latitude_min, extent$latitude_max),
add = TRUE)
# finish saving figure
dev.off()