R/aggregate_catchments.R

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

#' Aggregate Catchments
#' @description Aggregates catchments according to a set of outlet catchments. 
#' Network aggregation is completed using: See \code{\link{aggregate_network}}.
#' @param flowpath sf data.frame Flowpaths as generated by `refactor_nhdplus`
#' @param divide sf data.frame Reconciled catchment divides as generated by
#' `reconcile_catchment_divides`
#' @param outlets data.frame with "ID" and "type" columns. "ID" must be identifiers from
#' flowpath and divide data.frames. "type" should be "outlet", or "terminal".
#' "outlet" will include the specified ID.
#' "terminal" will be treated as a terminal node with nothing downstream.
#' @param zero_order list of vectors containing IDs to be aggregated into 0-order catchments.
#' @param coastal_cats sf data.frame with coastal catchments to be used with zero order.
#' @param da_thresh numeric See \code{\link{aggregate_network}}
#' @param only_larger boolean See \code{\link{aggregate_network}}
#' @param post_mortem_file rda file to dump environment to in case of error
#' @param keep logical passed along to \code{\link{clean_geometry}}

#' @details See \code{\link{aggregate_network}}
#'
#' @export
#' @importFrom sf st_is_empty st_drop_geometry st_as_sf
#' @importFrom dplyr filter left_join bind_rows select distinct
#' @examples
#' source(system.file("extdata", "walker_data.R", package = "hyRefactor"))
#' 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 = walker_fline_rec, 
#'                                    divide = walker_catchment_rec, 
#'                                    outlets = outlets)
#' plot(aggregated$cat_sets$geom, lwd = 3, border = "red")
#' plot(walker_catchment_rec$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)

aggregate_catchments <- function(flowpath, divide, outlets, zero_order = NULL,
                                 coastal_cats = NULL,
                                 da_thresh = NA, only_larger = FALSE, 
                                 post_mortem_file = NA, keep = NULL) {
  
  in_crs <- sf::st_crs(divide)
  
  if(any(grepl("MULTIPOLYGON", sf::st_geometry_type(divide, by_geometry = TRUE)))) {
    warning("found MULTIPOLYGONs in divides. They must be of type POLYGON.")
    divide <- clean_geometry(divide, ID = "ID", 
                             crs = in_crs, keep = NULL)
  }
  
  if(!is.null(zero_order)) {
    
    if(is.null(coastal_cats)) stop("must supply coastal_cats with zero order")
    
    if(st_crs(coastal_cats) != in_crs) st_transform(coastal_cats, in_crs)
    
    zero_flowpath <- filter(flowpath, member_COMID %in% do.call(c, zero_order))
  
    flowpath <- filter(flowpath, !ID %in% zero_flowpath$ID)
    
    coastal <- lapply(zero_order, function(x, coastal_cats, divide) {
      
      zero_cats <- filter(coastal_cats, FEATUREID %in% x)
      
      zero_div <- filter(divide, member_COMID %in% as.character(x))
      
      st_union(c(st_geometry(zero_cats), st_geometry(zero_div)))[[1]]
    }, coastal_cats = coastal_cats, divide = divide)
    
    coastal <- st_sfc(coastal, crs = in_crs)
    
    coastal <- st_sf(ID = names(coastal), geom = coastal)
  } else {
    coastal <- NULL
  }
  
  if(any(remove_head_div <- !flowpath$ID %in% flowpath$toID & 
         !flowpath$ID %in% divide$ID)) {
    remove_fpaths <- filter(flowpath, remove_head_div)
    message(paste("removing", nrow(remove_fpaths), 
                  "headwater/diversion flowlines without catchments."))
    flowpath <- filter(flowpath, !ID %in% remove_fpaths$ID)
  }

  agg_network <- aggregate_network(flowpath, outlets, da_thresh, only_larger, post_mortem_file)
  
  agg_network$cat_sets <- 
    unnest_flines(agg_network$cat_sets, 'set') %>% 
    left_join(select(divide, ID), by = c("set" = "ID")) %>% 
    st_as_sf() %>% 
    filter(!sf::st_is_empty(.)) %>%
    union_polygons(ID = "ID") %>% 
    clean_geometry(ID = "ID", 
                   crs = in_crs, 
                   keep = keep) %>% 
    select(.data$ID) %>% 
    left_join(agg_network$cat_sets, by = "ID")
  
  return(c(agg_network, list(coastal_sets = coastal)))
}