mr/hr subset function

hr_demo.R

subset_hr_mr <- function(lon, lat, out_dir) {
  
  dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
  
  point <- sf::st_sfc(sf::st_point(c(lon, lat)), 
                      crs = 4326) 
  
  mr_gpkg <- file.path(out_dir, "mr.gpkg")
  hr_gpkg <- file.path(out_dir, "hr.gpkg")
  out_gpkg <- file.path(out_dir, "out.gpkg")
  out_gpkg_internal <- file.path(out_dir, "internal_pecos.gpkg")
  
  outlet <- nhdplusTools::discover_nhdplus_id(point)
  
  png(file.path(out_dir, "mr.png"))
  
  mr <- nhdplusTools::plot_nhdplus(list(outlet), gpkg = mr_gpkg, 
                                   nhdplus_data = mr_gpkg, 
                                   overwrite = FALSE, 
                                   flowline_only = FALSE)
  
  dev.off()
  
  hu04 <- substr(dplyr::filter(mr$flowline, Hydroseq == min(Hydroseq))$REACHCODE, 
                 1, 4)
  
  hr_data <- nhdplusTools::download_nhdplushr(out_dir, hu04)
  
  hr <- nhdplusTools::get_nhdplushr(hr_data, out_gpkg = hr_gpkg, 
                                    overwrite = FALSE, check_terminals = TRUE)
  
  hr_outlet <- sf::st_intersection(hr$NHDPlusCatchment, 
                                   sf::st_transform(point, sf::st_crs(hr$NHDPlusCatchment)))
  hr_outlet_line <- dplyr::filter(hr$NHDFlowline, COMID == hr_outlet$FEATUREID)
  
  network <- nhdplusTools::get_UT(hr$NHDFlowline, hr_outlet$FEATUREID)
  
  subset <- nhdplusTools::subset_nhdplus(network, output_file = out_gpkg, 
                                         nhdplus_data = hr_gpkg, 
                                         return_data = TRUE, overwrite = TRUE)
  
  basin_boundary <- sf::st_sfc(sf::st_polygon(
    list(sf::st_geometry(sf::st_cast(mr$basin, "POLYGON"))[[1]][[1]])), 
    crs = sf::st_crs(mr$basin))
  
  outlets <- sf::st_intersection(dplyr::filter(hr$NHDFlowline, TerminalFl == 1), 
                                 sf::st_transform(basin_boundary, sf::st_crs(hr$NHDFlowline)))
  
  png(file.path(out_dir, "hr_outlets.png"))
  plot(sf::st_transform(basin_boundary, sf::st_crs(hr$NHDFlowline)))
  plot(sf::st_geometry(dplyr::filter(subset$NHDFlowline, TerminalFl == 1)), 
       add = TRUE, col = "red", lwd = 2)
  plot(sf::st_geometry(outlets), add = TRUE, col = "blue")
  dev.off()
  
  all_networks <- dplyr::filter(hr$NHDFlowline, TerminalPa %in% outlets$TerminalPa)$COMID
  
  internal <- nhdplusTools::subset_nhdplus(all_networks, 
                                           output_file = out_gpkg_internal, 
                                           nhdplus_data = hr_gpkg, 
                                           return_data = TRUE, overwrite = TRUE)

  png(file.path(out_dir, "hr.png"))

  plot(sf::st_simplify(sf::st_geometry(dplyr::filter(subset$NHDFlowline, StreamOrde > 2)), 
                       dTolerance = 0.005), col = "blue",
       lwd = (dplyr::filter(subset$NHDFlowline, StreamOrde > 2)$StreamOrde / 3))

  plot(sf::st_simplify(sf::st_geometry(dplyr::filter(internal$NHDFlowline, StreamOrde > 2)), 
                       dTolerance = 0.005), col = "red",
       lwd = (dplyr::filter(internal$NHDFlowline, StreamOrde > 2)$StreamOrde / 3), add = TRUE)
  
  dev.off()
  
}

Built with: https://usgs-r.github.io/nhdplusTools/index.html

Generates subset data and the following pngs:

NHDPlusV2.1 (MR) network.

Flowlines labeled as "terminal" in the HR subset.

HR network that connects to the outlet in blue -- does not connect in red.