Calculate level of nestedness of stream gages.

[dblodgett-usgs]

Some applications need to know if a stream gage is upstream of other stream gages for the purposes of working through the network in order of nesting level. e.g. calibrate parts of a model for headwaters only, then the regions just downstream of headwaters, then the next down, and so on. To accomplish this, an ordering of stream gages indicating their level of nesting is needed.

This level of nesting is actually very similar to stream order. Stream gages with nothing upstream get 1, stream gages with only 1s upstream get 2, stream gages with 2s upstream get 3, and so on.

The difference here is that, inlike stream order, "gage order" needs to be one more than the greatest upstream neighbor. This will require a tweak to the stream order logic but is much the same architecturally.

[dblodgett-usgs]

Draft implementation:

library(dplyr)
library(future)

get_dnid <- function(i, dnid) {
  
  out <- dnid[i]
  
  if(out == 0) {
    return()
  }
  
  c(out, get_dnid(out, dnid))
  
}

get_dn <- function(starts, to) {
  lapply(starts, \(x, to) {
    dn <- get_dnid(x, to)
    as.integer(tail(unlist(dn), 1))
  }, to)
}

network <- nhdplusTools::get_vaa(updated_network = TRUE)

network <- select(network, id = comid, toid = tocomid)

gages <- sf::read_sf("../geoconnex/ref_gages/out/usgs_gages.gpkg")

gages <- filter(gages, !is.na(nhdpv2_COMID) & nhdpv2_COMID %in% network$id)

gages <- left_join(gages, network, by = c("nhdpv2_COMID" = "id"))

gages <- select(gages, id = id, network_id = nhdpv2_COMID, network_toid = toid)

infer_connections <- function(x, network) {
  
  # these don't flow to anything -- order 1
  zero_ds <- x[x$network_toid == 0,]
  
  # these are only network hop and need to be handled out of band
  one_ds <- x[!x$id %in% zero_ds$id &
                            x$network_toid %in% x$network_id,]
  
  # we need to walk the network for all these
  to_trace <- x[!x$id %in% zero_ds$id &
                            !x$id %in% one_ds$id,]
  
  ### break the network downstream of gages we want to trace ###
  network$toid[network$id %in% to_trace$network_id] <- 0
  
  # get network ids
  hy_net <- hydroloom::make_index_ids(network)
  
  to_trace$indid <- hy_net$to_list$indid[match(to_trace$network_id, hy_net$to_list$id)]
  to_trace$to_indid <- hy_net$to_list$indid[match(to_trace$network_toid, hy_net$to_list$id)]
  
  all <- get_dn(to_trace$to_indid, hy_net$to)
  
  all[lengths(all) == 0] <- 0
  
  all <- as.integer(all)
  
  to_trace$down_id <- to_trace$id[match(all, to_trace$indid)]
  
  one_ds$down_id <- x$id[match(one_ds$network_toid, 
                                                x$network_id)]
  
  zero_ds$down_id <- rep("", nrow(zero_ds))
  
  out <- dplyr::bind_rows(zero_ds, one_ds, to_trace)
  
  out <- select(out, id = id, toid = down_id)
  
  stopifnot(all(x$id %in% out$id))
  
  out
}
  
net <- infer_connections(gages, network)

net <- sf::st_drop_geometry(net)

net <- hydroloom::sort_network(net)

index_ids <- hydroloom::make_index_ids(net)

froms <- hydroloom::make_fromids(index_ids)

order <- rep(1, length(froms$lengths))

for(i in seq_len(length(froms$lengths))) {
  
  l <- froms$lengths[i]
  
  # nothing to do if nothing upstream
  if(l > 0) {
    
    # these are the upstream orders
    orders <- order[froms$froms[1:l, i]]
    
    # Need the max upstream order for this work
    max_order <- max(orders)
    
    order[i] <- max_order + 1
    
  }
  
}

net$order <- order