--- title: "NHD data expedition day" author: "Nick Bruns" date: "2/18/2020" output: html_document --- Consult [these docs from this Download link](ftp://ftp.horizon-systems.com/NHDplus/NHDPlusV21/Documentation/NHDPlusV2_User_Guide.pdf) to understand NHD column values, their units etc. This is long, so just control-f on the information you need. [This map has all USGS gages](https://maps.waterdata.usgs.gov/mapper/index.html) [This tool is nice for visualizing HUCs](http://waterqualityexplorer.rc.duke.edu:3838/explorer/) ```{r setup, include=FALSE} # install.packages("nhdplusTools") #uncomment and run! library(nhdplusTools) library(tidyverse) library(sf) library(tmap) tmap_mode("view") #makes map plots interactive ``` --- ```{r target_hucs} #some gage ID's New_Hope_Creek_gage_id <- "02097314" Eno_gage_id <- "02085070" Yampa_gage_id <- "09244490" Yampa_gage_id_upstream <- "09239500" Weber_river_id <- "10128500" #UTah Salmon_river_ADK_id <- "04249200" #Aderondacks (good clean Hydrograph!) Salmon_river_ID_head_waters<- "13310700" Yellowstone_river_id <- "06214500" # at Billings Flint_river <- "02349605" #Georgia ``` ```{r code_demonstration} #First, we'll check how many machines work. So first, run this whole block, then we'll work through line by line. cur_gage_id <- Eno_gage_id cur_id_string<- paste0("USGS-",cur_gage_id) #just add prefix to string cur_site <- list(featureSource = "nwissite", #prepping below web quirie below, this just indicates the kind of data we have featureID = cur_id_string) cur_line <- navigate_nldi(cur_site, "UT", "") #returns all upstream flowlines cur_site_sf <- navigate_nldi(cur_site, "UT", "nwissite") #returns all upstream USGS gages #now map flowlines and gages tm_shape(cur_line) + tm_lines() + tm_shape(cur_site_sf) +tm_dots(col="red") ``` ```{r download_NHD_info} #Flowlines above are just the ID's . Download everything! save_file_name <- "Eno_River_NHD_subset.gpkg" #gkpg is an open source geospatial format if(!exists(save_file_name)){ subset_nhdplus(cur_line$nhdplus_comid, save_file_name, "download") } cur_nhd <- st_read(save_file_name) glimpse(cur_nhd) ``` ```{r make_nhd_plot} cur_density <- sum(cur_nhd$lengthkm)/max(cur_nhd$totdasqkm) #network density metric print(cur_density) #plot elevation cur_nhd %>% # filter(pathlength>2100) %>% #Necessary for Yampa river ggplot() + geom_line(aes(x=pathlength,y=minelevsmo,col=as.factor(streamorde))) + xlab("km to outlet") + ylab("elevation (cm)") + theme_dark() + scale_color_viridis_d("stream order") ``` ```{r wrap_all_above_code_into_functions} #great. That's it. Now we'll wrap all in functions for easier use. plot_basin <- function(gage_id_string){ id_string=paste0("USGS-",gage_id_string) site <- list(featureSource = "nwissite", featureID = id_string) print("getting site") discover_nldi_navigation(site) print("getting flowline") line <- navigate_nldi(site, "UT", "") #returns all upstream comids print("getting usptream sites") site_sf <- navigate_nldi(site, "UT", "nwissite") #returns all upstream gages! cur_plot <- tm_shape(line) + tm_lines() + tm_shape(site_sf) +tm_dots(col="red") print(cur_plot) return(line) } download_slim_nhd <- function(flow_line_with_comid,save_string){ subset_nhdplus(flow_line_with_comid$nhdplus_comid, save_string, "download") } plot_elevation <- function(cur_nhd,river_name){ cur_nhd %>% ggplot() + geom_line(aes(x=pathlength,y=minelevsmo,col=as.factor(streamorde))) + xlab("km to outlet") + ylab("elevation (cm)") + theme_dark() + scale_color_viridis_d("stream order") + ggtitle(river_name) } #Geomorphic Instantenous Unit Hydrograph plot_GIUH <- function(cur_nhd){ ggplot(cur_nhd) + geom_point(aes(y=totdasqkm,x=pathlength,col=as.character(streamorde))) + scale_color_viridis_d("stream order") + theme_dark() } compute_network_density<- function(cur_nhd){ sum(cur_nhd$lengthkm)/max(cur_nhd$totdasqkm) } plot_stream_order <- function(cur_nhd,matching_pal=T) { if(matching_pal){ cur_nhd %>% mutate(streamorde=as.factor(streamorde)) %>% tm_shape() + tm_lines(col="streamorde",lwd = 3,palette = "viridis") } else{ cur_nhd %>% mutate(streamorde=as.factor(streamorde)) %>% tm_shape() + tm_lines(col="streamorde",lwd = 3) } } ``` ```{r show_function_use_example} #setup: to try a different basin, # insert new code chunk (cmd + option + i) # paste this code chunk # modify this top part with a new name and USGS gage id. cur_river_name <- "New_Hope_Creek" cur_gage_id <- New_Hope_Creek_gage_id #or directly put a quoted gage id here print(cur_gage_id) #run code cur_file_name <- paste0(cur_river_name,"_NHD_subset.gpkg") cur_flow_line <- plot_basin(cur_gage_id) # download_slim_nhd(cur_flow_line,cur_file_name) #will fail if the file already exists cur_nhd <- st_read(cur_file_name) #load that NHD compute_network_density(cur_nhd) plot_elevation(cur_nhd,cur_river_name) plot_stream_order(cur_nhd ) plot_GIUH(cur_nhd ) ```