R/refactor_nhdplus.R

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

#' @title Refactor NHDPlus
#' @description A complete network refactor workflow has been packaged
#' into this function. Builds a set of normalized catchment-flowpaths from
#' input flowline features. See details and vignettes for more information.
#' @param nhdplus_flines data.frame raw nhdplus flowline features as
#' derived from the national seamless geodatabase.
#' @param split_flines_meters numeric the maximum length flowpath desired
#' in the output.
#' @param split_flines_cores numeric the number of processing cores to use
#' while splitting flowlines.
#' @param collapse_flines_meters numeric the minimum length of
#' inter-confluence flowpath desired in the output.
#' @param collapse_flines_main_meters numeric the minimum length of
#' between-confluence flowpaths.
#' @param out_refactored character where to write a geopackage containing
#' the split and collapsed flowlines.
#' @param out_reconciled character where to write a geopackage containing
#' the reconciled flowpaths.
#' @param three_pass boolean whether to perform a three pass collapse or
#' single pass.
#' @param purge_non_dendritic boolean passed on to prepare_nhdplus
#' @param exclude_cats integer vector of COMIDs to be excluded from collapse modifications.
#' @param events data.frame containing events as generated by nhdplusTools::get_flowline_index() 
#' @param warn boolean controls whether warning an status messages are printed
#' @details This is a convenient wrapper function that implements three phases
#' of the network refactor workflow: split, collapse, reconcile. See the
#' NHDPlus Refactor vignette for details of these three steps by running:
#' \code{vignette("refactor_nhdplus", package = "hyRefactor")}
#' @seealso
#' In addition to `prepare_nhdplus` from the nhdplusTools package,
#' The following three functions are used in the `refactor_nhdplus` workflow.
#' \enumerate{
#'   \item \code{\link{split_flowlines}}
#'   \item \code{\link{collapse_flowlines}}
#'   \item \code{\link{reconcile_collapsed_flowlines}}
#' }
#' @export
#' @importFrom nhdplusTools prepare_nhdplus
#' @examples
#' source(system.file("extdata",
#'                    "sample_flines.R",
#'                    package = "nhdplusTools"))
#' 
#' nhdplus_flowlines <- sf::st_zm(sample_flines)
#'
#' refactor_nhdplus(nhdplus_flines = nhdplus_flowlines,
#'                  split_flines_meters = 2000,
#'                  split_flines_cores = 2,
#'                  collapse_flines_meters = 500,
#'                  collapse_flines_main_meters = 500,
#'                  out_refactored = "temp.gpkg",
#'                  out_reconciled = "temp_rec.gpkg",
#'                  three_pass = TRUE,
#'                  purge_non_dendritic = FALSE,
#'                  warn = FALSE)
#'                  
#' unlink("temp.gpkg")
#' unlink("temp_rec.gpkg")

refactor_nhdplus <- function(nhdplus_flines,
                             split_flines_meters,
                             split_flines_cores,
                             collapse_flines_meters,
                             collapse_flines_main_meters,
                             out_refactored,
                             out_reconciled,
                             three_pass = FALSE,
                             purge_non_dendritic = TRUE,
                             exclude_cats = NULL,
                             events = NULL,
                             warn = TRUE) {

  if ("FTYPE" %in% names(nhdplus_flines)) {
    nhdplus_flines <- dplyr::inner_join(
      select(sf::st_zm(nhdplus_flines), COMID, REACHCODE, FromMeas, ToMeas),
      sf::st_drop_geometry(nhdplus_flines) %>%
        nhdplusTools::prepare_nhdplus(0, 0, 0, purge_non_dendritic = purge_non_dendritic,
                        warn = warn), by = "COMID") %>%
      sf::st_as_sf()
  }

  in_proj <- sf::st_crs(nhdplus_flines)

  flines <- nhdplus_flines %>%
    sf::st_cast("LINESTRING", warn = warn) %>%
    sf::st_transform(5070) %>%
    split_flowlines(split_flines_meters, para = split_flines_cores, 
                    avoid = exclude_cats, events = events)

  rm(nhdplus_flines)

  if (warn) {
    message("flowlines split complete, collapsing")
  }
  
  exclude_cats <- c(exclude_cats, dplyr::filter(flines, !is.na(event_REACH_meas))$COMID,
                    dplyr::filter(flines, !is.na(event_REACH_meas))$toCOMID)

  if (three_pass) {
    collapsed_flines <-
      collapse_flowlines(sf::st_set_geometry(flines, NULL),
                         (0.25 * collapse_flines_meters / 1000),
                         TRUE,
                         (0.25 * collapse_flines_main_meters / 1000),
                         exclude_cats)

    collapsed_flines <-
      collapse_flowlines(collapsed_flines,
                         (0.5 * collapse_flines_meters / 1000),
                         TRUE,
                         (0.5 * collapse_flines_main_meters / 1000),
                         exclude_cats,
                         warn = FALSE)

    collapsed_flines <-
      collapse_flowlines(collapsed_flines,
                         (collapse_flines_meters / 1000),
                         TRUE,
                         (collapse_flines_main_meters / 1000),
                         exclude_cats,
                         warn = FALSE)
  } else {
    collapsed_flines <-
      collapse_flowlines(sf::st_set_geometry(flines, NULL),
                         (collapse_flines_meters / 1000),
                         TRUE,
                         (collapse_flines_main_meters / 1000),
                         exclude_cats)
  }

  select(flines, COMID) %>%
    dplyr::inner_join(collapsed_flines, by = "COMID") %>%
    sf::st_as_sf() %>%
    sf::st_transform(in_proj) %>%
    sf::st_write(out_refactored, layer_options = "OVERWRITE=YES",
                 quiet = !warn)

  if (warn) {
    message("collapse complete, out collapse written to disk, reconciling")
  }

  collapsed <- reconcile_collapsed_flowlines(collapsed_flines,
                                             select(flines, COMID),
                                             id = "COMID")

  collapsed[["member_COMID"]] <-
    unlist(lapply(collapsed$member_COMID,
                  function(x) paste(x, collapse = ",")))

  sf::st_write(sf::st_transform(collapsed, in_proj),
               out_reconciled,
               layer_options = "OVERWRITE=YES",
               quiet = !warn)
}