initial 02_process_data.R and 03_make_figures.R

01_download_data.R

@@ -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)

02_process_data.R

@@ -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)
+}

03_make_figures.R

@@ -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()