NHDPlusHR with hygeo output

dblodgett@usgs.gov

First, we’ll load up three data sets. 1. the mainstem identifiers from NHDPlusV2 (LevelPathIDs) 2. NHDPlusHR Flowlines which have the NHDPlusV2 LevelPathID as an attribute. 3. Reconciled flowline output from the processing shown in vignette("sugar_creek_refactor")

library(nhdplusTools)
#> USGS Support Package: https://owi.usgs.gov/R/packages.html#support
library(hyRefactor)
#> USGS Support Package: https://owi.usgs.gov/R/packages.html#support
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidyr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(hygeo)
src_gpkg <- system.file("gpkg/sugar_creek_fort_mill.gpkg", package = "hygeo")

mr_lp <- st_drop_geometry(read_sf(src_gpkg, "NHDFlowline_Network")) %>%
  select(COMID, LevelPathI)

hr_fline <- read_sf(src_gpkg, "hr_NHDPlusFlowline") %>%
  select(NHDPlusID = COMID, LevelPathI, Hydroseq, mr_LevelPathI)

fline <- read_sf("../tests/testthat/data/sugar_creek_hyRefactor.gpkg", "reconcile")

waterbody_edge_list <- get_waterbody_edge_list(fline,
                                               waterbody_prefix = "fp-")
waterbody <- get_flowpath_data(fline,
                                waterbody_edge_list,
                                catchment_prefix = "cat-")

Given these three datasets we need to do a little manipulation to get the identifiers straight.

First we create a cross walk for mainstem ids. This is a little awkward. See: https://github.com/dblodgett-usgs/hyRefactor/issues/5 Using the cross walk, we can add attributes generate outlet node locations from the NHDPlusHR flowlines and add the mainstem identifier cross walk to them.

main_id_xwalk <- select(st_drop_geometry(waterbody),
                        local_id = ID, main_id) %>%
  left_join(get_nhd_crosswalk(fline, catchment_prefix = "cat-"), by = "local_id") %>%
  mutate(COMID = as.integer(COMID)) %>%
  select(-local_id) %>%
  distinct() %>%
  left_join(mr_lp, by = "COMID") %>%
  select(-COMID)

hr_nodes <- st_sf(st_drop_geometry(hr_fline),
                  geom =  st_geometry(nhdplusTools::get_node(hr_fline))) %>%
  left_join(main_id_xwalk, by = c("mr_LevelPathI" = "LevelPathI"))

Now we just need to make sure things are in the project we want to use for analysis, the columns of the input points are right, and pass them in.

hyl <- st_transform(select(hr_nodes, NHDPlusID, main_id), 5070)

waterbody <- st_transform(waterbody, 5070)

(hyl_out <- get_hydrologic_locaton(hyl, waterbody))
#> Simple feature collection with 5053 features and 4 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 1350590 ymin: 1429382 xmax: 1378126 ymax: 1468134
#> projected CRS:  NAD83 / Conus Albers
#> First 10 features:
#>         NHDPlusID    ID  measure    offset                    geom
#> 1  15001500020498 cat-4  7.94341  1.243441 POINT (1356654 1449829)
#> 2  15001500020498 cat-4  7.94341  1.243441 POINT (1356654 1449829)
#> 3  15001500058426 cat-4 31.43707 18.909778 POINT (1355922 1451139)
#> 4  15001500058426 cat-4 31.43707 18.909778 POINT (1355922 1451139)
#> 5  15001500058428 cat-4 92.70761 52.079451 POINT (1354513 1454544)
#> 6  15001500058428 cat-4 92.70761 52.079451 POINT (1354513 1454544)
#> 7  15001500058429 cat-4 60.03731 16.278611 POINT (1355446 1452818)
#> 8  15001500058429 cat-4 60.03731 16.278611 POINT (1355446 1452818)
#> 9  15001500135210 cat-4 48.21983 24.949110 POINT (1355369 1452009)
#> 10 15001500135210 cat-4 48.21983 24.949110 POINT (1355369 1452009)

Now we can map up the matches we got.

hr_matched <- filter(hr_fline, NHDPlusID %in% hyl_out$NHDPlusID)

mapview::mapview(fline, lwd = 4, color = "blue", layer.name = "NHDPlusV2") + 
  mapview::mapview(hr_matched, lwd = 3, color = "red", layer.name = "NHDPlusHR") +
  mapview::mapview(hyl_out, cex = 2, color = "black", layer.name = ("NHDPlusHR Outlets"))