From 7eddee51b4ea1f1c9512ab9169981468801cabb3 Mon Sep 17 00:00:00 2001 From: Ben Varick Date: Mon, 16 Oct 2023 12:17:24 -0500 Subject: [PATCH] initial 02_process_data.R and 03_make_figures.R --- 01_download_data.R | 12 +++++++++--- 02_process_data.R | 16 ++++++++++++++++ 03_make_figures.R | 47 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 72 insertions(+), 3 deletions(-) diff --git a/01_download_data.R b/01_download_data.R index 2b8eb61..1e2fc02 100644 --- a/01_download_data.R +++ b/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) \ No newline at end of file diff --git a/02_process_data.R b/02_process_data.R index e69de29..8c877cc 100644 --- a/02_process_data.R +++ b/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) +} diff --git a/03_make_figures.R b/03_make_figures.R index e69de29..25813be 100644 --- a/03_make_figures.R +++ b/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()