local <- (Sys.getenv("BUILD_VIGNETTES") == "TRUE") knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=6, fig.height=4, eval=local ) options(rmarkdown.html_vignette.check_title = FALSE) options(scipen = 9999) options("rgdal_show_exportToProj4_warnings"="none")
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)
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)
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)
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:
There are two caveats here.
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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.