Experiment to find main-flowpaths through catchment aggregates.

demo_aggregation.R

# start at upstream most head water nexus and trace downstream
# find groups bound by outlets and 
# mark visited catchments in catchment topology vector
# keep going until all outlets and headwaters have been consumed
demo_aggregation <- function() {
library(dplyr)
library(sf)

load("testing.rda")

outlet_points <- filter(flowpath, ID %in% outlets$ID) %>%
  nhdplusTools::get_node()

# Make sure flowpath is sorted correctly
flowpath <- arrange(flowpath, desc(Hydroseq))

# Do some id redirection for performance.
orig_ids <- select(sf::st_drop_geometry(flowpath), ID, toID)
orig_ids["id"] <- seq(1, nrow(flowpath))

flowpath$toid <- match(flowpath$toID, flowpath$ID)
flowpath$id <- seq(1, nrow(flowpath))

outlets <- left_join(outlets, select(orig_ids, -toID), by = "ID")
outlets$set <- 1:nrow(outlets)

# We'll start at all these
heads <- flowpath$id[!flowpath$id %in% flowpath$toid]

# We can stop once these conditions have been met.
outlet_count <- nrow(outlets)
headwater_count <- length(heads)

# need a little network walker function.
get_dwn <- function(ID, toid) {
  next_dn <- toid[ID]
  if(is.na(next_dn)) {
    return(ID)
  } else {
    return(c(ID, get_dwn(next_dn, toid)))
  }
}

# Set up counters for the while loop
o_c <- 1
h_c <- 1

while(o_c < outlet_count & h_c < headwater_count) {
  head <- heads[h_c]
  
  path <- get_dwn(head, flowpath$toid)
  
  sets <- outlets$id %in% path
  
  breaks <- outlets[sets, ]
  nr_breaks <- nrow(breaks)
  
  plot(sf::st_geometry(flowpath))
  plot(sf::st_geometry(flowpath[is.na(flowpath$toid), ]), add = TRUE, lwd = 3)
  plot(sf::st_geometry(outlet_points), add = TRUE)
  plot(sf::st_geometry(flowpath[flowpath$id == head, ]), add = TRUE, col = "blue", lwd = 2)

  # If this path only goes to one outlet, it is within an aggregate
  ## do nothing. Or maybe we can use this for aggregation later?
  # If this path goes through more than one outlet, it contains flowpaths.
  if(nr_breaks > 1) {
    paths <- split(path, 
                   cut(path, 
                       breaks = c(0, breaks$id), 
                       labels = c(breaks$set)))
    
    # the top one isn't useful.
    paths <- paths[2:length(paths)]
    
    fline_sets$set[as.integer(names(paths))] <- 
      lapply(paths, function(x) flowpath$ID[x])
    
    flowpath$toid[unlist(paths)] <- NA
    
    o_c <- o_c + (nr_breaks - 1)
    
    # in the current path removal
    plot(sf::st_geometry(flowpath[unlist(paths), ]), add = TRUE, col = "red")
    
    
  } 
  
  h_c <- h_c + 1
  
}

plot(sf::st_geometry(flowpath))
plot(sf::st_geometry(flowpath[is.na(flowpath$toid), ]), add = TRUE, lwd = 3)
plot(sf::st_geometry(outlet_points), add = TRUE)

set_length <- lengths(fline_sets$set)

head_outlets <- outlets[set_length == 0, ] %>%
  left_join(lps, by = "ID") 

plot(sf::st_geometry(flowpath[flowpath$ID %in% head_outlets$ID, ]), lwd = 2, col = "red", add = TRUE)

head_outlets <- head_outlets %>%
  select(LevelPathID, head_out_Hydroseq = Hydroseq)

head_paths <- filter(flowpath, LevelPathID %in% head_outlets$LevelPathID) %>%
  left_join(head_outlets, by = "LevelPathID") %>%
  filter(Hydroseq >= head_out_Hydroseq)
  
plot(sf::st_geometry(head_paths), col = "blue", add = TRUE)
}

gifski::save_gif(demo_aggregation())

Black is the core mainstem network. Red are first order outlets. Blue are mainstem paths for first order aggregates.