Added more comments and specified library::function calls

NHD_navigate.R

@@ -8,8 +8,7 @@
 ##         in this example
 ## NEEDED MODS
 #  0 - a) create setup script for downloading and formatting data 
-#      b) downloading data from SB into the proper formats
-#      c) ensuring the proper libraries are installed and loaded.
+#      b) ensuring the proper libraries are installed and loaded.
 #  1 - Functionalize script with input arguments to cover workspace, results, etc.
 #  2 - Add Gages to workflow
 #  3 - Add NID/Waterbodies COMIDs to workflow
@@ -17,63 +16,79 @@
 #  4 - Spit out segment length, travel time
 #  5 - identify incremental catchments associated with each downstream POI
 ##################################################
+
 # Network navigation for upstream/downstream from a COMID of interest
 NetworkNav<-function(inCom,nhdDF,navType){
   if (navType=="up"){
-    Seg<-get_UM(nhdDF,inCom,include=T)
+    # Upstream network navigation
+    Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
   }else{
-    Seg<-get_DM(nhdDF,inCom,distance=1050,include=T)
+    # Downstream network navigation (distance currently hard-coded)
+    Seg<nhdplusTools::get_DM(nhdDF,inCom,distance=1050,include=T)
   }
   return(Seg)
 }
 
-library(nhdplusTools)
-library(tidyverse)
-library(sf)
-library(rgdal)
+# POI Creation
+POI_creation<-function(inSegs){
+  # break segs it to points
+  Segpts <- st_geometry(inSegs) %>% 
+    lapply(., function(x) {sf::st_sfc(x) %>% sf::st_cast(., 'POINT')})
+  # subset the last point from each geometry, make a POINT sf object
+  Seg_ends <- sapply(Segpts, function(p) {
+    p[c(length(p))]}) %>% sf::st_sfc() %>% sf::st_sf('geom' = .)
+  return(Seg_ends)
+}
+
+library(nhdplusTools) #https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
+library(tidyverse) # a collection of libraries that work with eacher other (dplyr)
+library(sf) # simple features - a library for handling spatial data
+library(rgdal) # rgdal - an open source GIS library that's more akin to ESRI and handles shapefile
 
-#https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
+# Set working directory
 setwd("D:/abock/NHM_FY20/gfv2.0")
 
 #######################
-# Read Region 06 HUC12s
-HUC12<-readOGR("data/hu_points.gpkg","hu_points")
+# Read Region 06 HUC12 outlets, get this data at https://www.sciencebase.gov/catalog/item/5dbc53d4e4b06957974eddae
+HUC12<-rgdal::readOGR("data/hu_points.gpkg","hu_points")
+# Add new region field
 HUC12@data$Region<-substr(HUC12@data$HUC12,1,2)
-comIDs<-HUC12@data %>% filter(Region == "06") %>% select (COMID)
+# Get the COMIDS of HUC12 outlets within Region 06
+comIDs<-HUC12@data %>% dplyr::filter(Region == "06") %>% dplyr::select (COMID)
 
 # Set path to national NHDPlus v2 seamless data; 
 Print ("This will take awhile, go take a walk")
 if (!file.exists("data/nhdplus_flowline_attributes.rds")){ 
-  download_nhdplusv2("data")
-  nhdplus_path("data/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
-  stage_national_data(include=c("attribute","flowline"),output_path="data") # Get catchments later?
+  # Download seamless data to this location
+  nhdplusTools::download_nhdplusv2("data")
+  # Set the path
+  nhdplusTools::nhdplus_path("data/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
+  # Stage the data as RDS for easy access
+  nhdplusTools::stage_national_data(include=c("attribute","flowline","catchment"),output_path="data") # Get catchments later?
 }else{
-  print ("NHDPlusv2 data downloaded and staged")
+  print ("NHDPlusv2 data already downloaded and staged")
 }
 
-# Read flowline attribute data frames
+# Read flowline attribute data frames for entire seamless database
 nhdDF<-readRDS("data/nhdplus_flowline_attributes.rds") # %>% select(COMID, Pathlength, LENGTHKM, LevelPathI, UpHydroseq, Hydroseq, DnHydroseq)
 # Read flowline simple features
 nhd<-readRDS("data/nhdplus_flowline.rds") 
 # Return flowlinest that are POIs
-huc12POIsegs<-nhd %>% filter(COMID %in% comIDs[,1]) # 15 missing?
+huc12POIsegs<-nhd %>% dplyr::filter(COMID %in% comIDs[,1]) # 15 missing?
 
 #########################
-# break segs it to points
-Segpts <- st_geometry(huc12POIsegs) %>% 
-  lapply(., function(x) {st_sfc(x) %>% st_cast(., 'POINT')})
-# subset the last point from each geometry, make a POINT sf object
-Seg_ends <- sapply(Segpts, function(p) {
-  p[c(length(p))]}) %>% st_sfc() %>%st_sf('geom' = .)
+# Create POIs from HUC12 segments
+huc12POIs<-POI_creation(huc12POIsegs)
 # Write out shapefile representing POIs for given theme
-st_write(Seg_ends,dsn="workspace/huc12POIs.shp",delete_dsn = T)
-write_rds(Seg_ends,path="workspace/huc12POIS.rds")
+sf::st_write(huc12POIs,dsn="workspace/huc12POIs.shp",delete_dsn = T)
+write_rds(huc12POIs,path="workspace/huc12POIS.rds")
 
 ##########################
-# Navigate upstream along the mainstems
+# Navigate upstream along the mainstems from HUC12 POIS
 upNet<-unique(unlist(lapply(comIDs,NetworkNav,nhdDF,"up")))
 # Subset by dangling segments that need downstream navigation to be connected
-upNet_DF<-nhd %>% filter(COMID %in% upNet) %>% 
+#        reduces the number of features we need to connect
+upNet_DF<-nhd %>% dplyr::filter(COMID %in% upNet) %>% 
   filter(!DnHydroseq %in% Hydroseq)
 # Nagivate downstream along the mainstems
 downNet<-unique(unlist(lapply(upNet_DF$COMID,NetworkNav,nhdDF,"down")))
@@ -82,20 +97,21 @@ huc12Segs<-unique(c(upNet,downNet))
 
 # Subset nhd segs to navigation results and write to shapefile
 huc12Segs<-nhd %>% filter(COMID %in% huc12Segs)
-st_write(huc12Segs,dsn="workspace/huc12segs.shp",delete_dsn = T)
+sf::st_write(huc12Segs,dsn="workspace/huc12segs.shp",delete_dsn = T)
 write_rds(Seg_ends,path="workspace/huc12segs.rds")
 
+##########################
 # Derive confluence POIs where they are absent
-nhdConf<-nhd %>% filter(COMID %in% huc12Segs$COMID)
-confluences<-nhdConf %>% group_by(DnHydroseq) %>% filter(n()>1) %>% filter(!COMID %in% huc12POIsegs$COMID)
+nhdConf<-nhd %>% dplyr::filter(COMID %in% huc12Segs$COMID) # check on redundancy here
+
+confluences<-nhdConf %>% dplyr::group_by(DnHydroseq) %>% dplyr::filter(n()>1) %>% 
+  dplyr::filter(!COMID %in% huc12POIsegs$COMID)
   
-confpts <- st_geometry(confluences) %>% 
-  lapply(., function(x) {st_sfc(x) %>% st_cast(., 'POINT')})
-# subset the last point from each geometry, make a POINT sf object
-conf_ends <- sapply(confpts, function(p) {
-  p[c(length(p))]}) %>% st_sfc() %>%st_sf('geom' = .)
+# Create POIs from HUC12 segments
+huc12POIs<-POI_creation(confluences)
+
 # Write out shapefile representing POIs for given theme
-st_write(conf_ends,dsn="workspace/confPOIs.shp",delete_dsn = T)
+sf::st_write(conf_ends,dsn="workspace/confPOIs.shp",delete_dsn = T)
 write_rds(conf_ends,path="workspace/confPOIs.rds")
 
 ##########################################