In dblodgett-usgs/hyRefactor: Hydrologic Network Refactoring Tools Based on HY_Features

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6,
  fig.height=4
)
options(rmarkdown.html_vignette.check_title = FALSE)
options(scipen = 9999)
options("rgdal_show_exportToProj4_warnings"="none")

Refactoring Catchments with NHDPlusTools

Load Network and Refactor

For this example, we will use some data available from the nhdplusTools package for testing and examples. fac_sample, fdr_sample, flowline_sample, catchment_sample are all included as sample data in nhdplusTools.

Before we can do anything with the catchments, we need to run the NHDPlus Refactor workflow documented over in vignette("refactor_nhdplus").

library(sf)
library(nhdplusTools)
library(hyRefactor)

source(system.file("extdata", "walker_data.R", package = "hyRefactor"))

fac <- walker_fac@ptr$filenames
fdr <- walker_fdr@ptr$filenames

ref <- tempfile(fileext = ".gpkg")
rec <- tempfile(fileext = ".gpkg")

refactor_nhdplus(nhdplus_flines = walker_flowline,
                 split_flines_meters = 2000,
                 collapse_flines_meters = 1000,
                 collapse_flines_main_meters = 1000,
                 split_flines_cores = 2,
                 out_refactored = ref,
                 out_reconciled = rec,
                 three_pass = TRUE,
                 purge_non_dendritic = FALSE,
                 warn = FALSE)

flowline_ref <- read_sf(ref)
flowline_rec <- read_sf(rec)

Let's start by looking at a single sample catchment. In the first plot, you can see the input is a single polygon and a set of catchment flowlines. For reference, the Flow Direction Raster (FDR) and Flow Accumulation Raster (FAC) are also plotted.

sample_catchment <- dplyr::filter(walker_catchment, FEATUREID == 5329435)

sample_flowline <- dplyr::filter(flowline_ref, as.integer(COMID) == 5329435)

plot(st_geometry(sample_catchment))
plot(sample_flowline["COMID"], add = TRUE)

plot(walker_fdr)
plot( walker_fac)

Now we can run the split_catchment_divide() function which is designed to take one catchment and its associated (split) flowlines. The split flowlines are the "refactored" but not "reconciled" output of the nhdplus_refactor() function we ran above. Here we run the split_catchment() function and plot up the resulting data to show what it did.

sample_catchment <- sf::st_transform(sample_catchment, 
                                     sf::st_crs(walker_fdr))
sample_flowline <- sf::st_transform(sample_flowline, 
                                     sf::st_crs(walker_fdr))

split_cat <- split_catchment_divide(sample_catchment, sample_flowline, walker_fdr,  walker_fac)

plot(st_geometry(split_cat), col = NA, border = "red")
plot(st_geometry(sample_catchment), add = TRUE)
plot(sample_flowline["COMID"], lwd = 2, add= TRUE)

As you can see, the flowline in question was split into five pieces by nhdplus_refactor() and the cooresponding catchment (black) was similarly broken up into 5 sub-catchments (red).

The split_catchment_divide() function can be run against a set of catchments using the reconcile_catchment_divides() function. This functiona can call split_catchment() in parallel, unions catchments according to the output of nhdplus_refactor(), and assembles the results back into an sf data.frame.

walker_catchment <- sf::st_transform(walker_catchment, sf::st_crs(walker_fdr))
flowline_ref <- sf::st_transform(flowline_ref, sf::st_crs(walker_fdr))
flowline_rec <- sf::st_transform(flowline_rec, sf::st_crs(walker_fdr))

split_cats <- reconcile_catchment_divides(catchment = walker_catchment, 
                                          fline_ref = flowline_ref, 
                                          fline_rec = flowline_rec, 
                                          fdr = fdr, 
                                          fac = fac, 
                                          para = 4)

plot(st_geometry(split_cats), col = NA, border = "red")
plot(st_geometry(walker_catchment), col = NA, border = "black")
plot(st_geometry(split_cats), col = NA, border = "red")
plot(st_geometry(walker_catchment), col = NA, border = "black", add = TRUE)

Aggregation

With our flowline network and catchments all refactored and the geometry reconciled, we could also aggregate the resulting network to a selected set of outlet locations.

outlets <- data.frame(ID = c(31, 3, 5, 1, 45, 92),
                      type = c("outlet", "outlet", "outlet", "terminal", "outlet", "outlet"),
                      stringsAsFactors = FALSE)

aggregated <- aggregate_catchments(flowpath = flowline_rec, 
                                   divide = split_cats,
                                   outlets = outlets)

plot(aggregated$cat_sets$geom, lwd = 3, border = "red")
plot(split_cats$geom, lwd = 1.5, border = "green", col = NA, add = TRUE)
plot(walker_catchment$geom, lwd = 1, add = TRUE)
plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE)
#'
plot(aggregated$cat_sets$geom, lwd = 3, border = "black")
plot(aggregated$fline_sets$geom, lwd = 3, col = "red", add = TRUE)
plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE)

Splitting Catchments

The split_catchment_divide() function can be used in two ways. If provided with multiple flowpaths per catchment boundary, it will split catchments along their length. If provided one or more with the lr parameter set to TRUE, it will also split each catchment into two pieces along the flowpath. In headwater catchments, the left-right split will extend upstream along the maximum upstream flow accumulation path till it reaches the catchment divide. split_catchment_divide() is called in the refactor_nhdplus() workflow but this left-right split is not applied then. This next block shows how to apply the split to aggregated results.

NOTE: left / right split is not function in v0.4.8

split_catchments <- do.call(c, lapply(c(1:nrow(aggregated$cat_sets)), 
                                      function(g, ac, af, fdr, fac) {
                                        split_catchment_divide(ac[g, ], af[g, ], 
                                                               fdr, fac, lr = FALSE)
                                      }, 
                                      ac = aggregated$cat_sets, 
                                      af = aggregated$fline_sets, 
                                      fdr = walker_fdr, 
                                      fac = walker_fac))

Here we see the split geometry in black and the original catchments in red.

plot(split_catchments, lwd = 3)
plot(st_geometry(aggregated$cat_sets), lwd = 2, border = "red", add = TRUE)

Here we see the split geometry in black again with the flowline used for the split in red.

plot(split_catchments, lwd = 3)
plot(aggregated$fline_sets$geom, lwd = 1, col = "red", add = TRUE)

Identifiers

So you want to know how the identifiers from the original NHDPlus Catchments relate to the output? Let's look at what's in the output.

(input_ids <- walker_flowline$COMID)

(refactored_ids <- flowline_rec$ID)

(refactored_id_mapping <- flowline_rec$member_COMID)

(aggregated_ids <- aggregated$cat_sets$ID)

(aggregated_id_mapping <- aggregated$cat_sets$set)

This may look like a complicated mess but there's structure here. Stepping through the identifiers, we have:

1. NHDPlusV2 COMIDs (shared between catchment divides and flowlines)
2. Refactored Catchment Identifiers (shared between catchment divides and flowpaths)

There are two caveats here.

1. Because the source NHDPlusV2 catchments were split, the relationship between refactored catchments and source catchments requires a sequence (.1, .2, .. .10, etc. -- upstream to downstream) to differentiate the parts.
2. Aggregated catchments use the identifer of the outlet catchment.

Given these caveats, we can build a complete lookup table from source catchment to output aggregate identifier.

refactor_lookup <- dplyr::select(st_drop_geometry(flowline_rec), ID, member_COMID) %>%
  dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
  hyRefactor:::unnest_flines(col = "member_COMID") %>%
  dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
  dplyr::rename(reconciled_ID = ID)

aggregate_lookup_fline <- dplyr::select(st_drop_geometry(aggregated$fline_sets), ID, set) %>%
  hyRefactor:::unnest_flines() %>%
  dplyr::rename(aggregated_flowline_ID = ID, reconciled_ID = set)

aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(aggregated$cat_sets), ID, set) %>%
  hyRefactor:::unnest_flines() %>%
  dplyr::rename(aggregated_catchment_ID = ID, reconciled_ID = set)

(lookup_table <- tibble::tibble(NHDPlusV2_COMID = input_ids) %>%
  dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID") %>%
  dplyr::left_join(aggregate_lookup_fline, by = "reconciled_ID") %>%
  dplyr::left_join(aggregate_lookup_catchment, by = "reconciled_ID"))

Ta Da!! Have fun and report bugs here.

dblodgett-usgs/hyRefactor documentation built on Aug. 25, 2023, 9:09 p.m.