Strahler stream order using igraph or sfnetwork in R

Question

I cannot fathom how to derive strahler order in R. Here's an example in postgres and neo4j. An attempt in R

There are three rules (from the GRASS 7.8 Manual):

• if the node has no children, it's Strahler order is 1.
• if the node has one and only one tributary with Strahler greatest order i, and all other tributaries have order less than i, then the order remains i.
• if the node has two or more tributaries with greatest order i, then the Strahler order of the node is i + 1.

Here's what I would expect

library(sfnetworks)
library(igraph)
library(sf)
library(dplyr)
library(tidygraph)
library(RColorBrewer)

# Create an example network.
n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(1, 2)))
n03 = st_sfc(st_point(c(1, 3)))
n04 = st_sfc(st_point(c(1, 4)))
n05 = st_sfc(st_point(c(2, 1)))
n06 = st_sfc(st_point(c(2, 3)))
n07 = st_sfc(st_point(c(2, 4)))
n08 = st_sfc(st_point(c(3, 2)))
n09 = st_sfc(st_point(c(3, 3)))
n10 = st_sfc(st_point(c(3, 4)))
n11 = st_sfc(st_point(c(4, 2)))
n12 = st_sfc(st_point(c(4, 4)))

from = c(1, 2, 2, 3, 3, 5, 5, 8, 8, 9, 9)
to = c(5, 3, 6, 4, 7, 2, 8, 9, 11, 10, 12)

nodes = st_as_sf(c(n01, n02, n03, n04, n05, n06, n07, n08, n09, n10, n11, n12))
edges = data.frame(from = from, to = to)

G = sfnetwork(nodes, edges) %>%
  convert(to_spatial_explicit, .clean = TRUE)

nodes = st_as_sf(G, "nodes")
edges = st_as_sf(G, "edges")

# expected order
edges$expected_order = c(4,2,1,1,1,3,3,2,1,1,1) 

cols = brewer.pal(4, "Blues")
pal = colorRampPalette(cols)

plot(st_geometry(edges))
plot(edges["expected_order"], 
     lwd = 4, , 
     add = TRUE,
     col = pal(4)[edges$expected_order])
legend(x = "topright",
       legend = c("4","3","2","1"),          
       lwd = 4,
       col = pal(4)[edges$expected_order],
       title = "strahler order")
plot(nodes, pch = 20, add = TRUE)

Here's what I tried curtesy of jsta/streamnet/stream_order.R, which I can't load due to missing packages

stream_order_igraph <- function(tree){
  
  tree <- as.igraph(tree)
  
  leaf_nodes <- which(degree(tree,
                                   v = igraph::V(tree),
                                   mode = "in") == 0,
                            useNames = TRUE)
  
  base_order <- 1
  
  edgelist   <- data.frame(as_edgelist(tree))
  edgelist$order <- NA
  names(edgelist)[c(1,2)] <- c("from", "to")
  edgelist$order[edgelist$from %in% leaf_nodes] <- base_order
  
  tree <- igraph::delete.vertices(tree, leaf_nodes)
  
  while(igraph::vcount(tree) >= 1){
    base_order <- max(edgelist$order, na.rm = TRUE) + 1
    leaf_nodes <- which(degree(tree, v = igraph::V(tree),
                                     mode = "in") == 0,
                              useNames = TRUE)
    
    raised_nodes <- sapply(leaf_nodes,
                           function(x) all(edgelist$order[edgelist$to == x] == base_order - 1))
    raised_nodes <- which(raised_nodes)
    flat_nodes <- leaf_nodes[!(leaf_nodes %in% raised_nodes)]
    
    edgelist$order[edgelist$from %in% raised_nodes] <- base_order
    edgelist$order[edgelist$from %in% flat_nodes] <- base_order - 1
    
    tree <- igraph::delete.vertices(tree, leaf_nodes)
    
  }
  edgelist$order
}

stream_order_igraph(G)

> stream_order_igraph(G)
 [1]  4  3  3  3  3  2  2 NA NA NA NA

Answer 1

I have found a solution that converts class igraph to class phylo then uses phytools::StrahlerNumber. I had to modify phytools::igraph_to_phylo and reverse the order of my edges to get it to work.

library(phytools)
library(igraph)
library(sfnetworks)
library(sf)
library(dplyr)
library(RColorBrewer)

# reverse the edge direction 
transposeGraph <- function(g) {
  g %>% get.edgelist %>%
    {cbind(.[, 2], .[, 1])} %>%
    graph.edgelist
}

# convert igraph class to phylo class
# and calculate strahler number 
igraphStrahler <- function(g){
  
  if (!igraph::is_simple(g) |
      !igraph::is_connected(g) |
      !igraph::is_dag(g)) {
    stop("Taxon graph is not a simple, connected, directed acylic graph")
  }
  
  root = which(sapply(V(g), 
                      function(x) length(neighbors(g, x, mode = "in"))) == 0)
  leaves = which(sapply(V(g), 
                          function(x) length(neighbors(g, x, mode = "out"))) == 0)
  
  g <- g %>% 
    set_vertex_attr("leaf", index = leaves, TRUE) %>%
    set_vertex_attr("root", index = root, TRUE)
  
  traverse <- igraph::dfs(g, root)
  
  is_leaf <- igraph::vertex_attr(g, "leaf", traverse$order)
  is_leaf[which(is.na(is_leaf))] <- FALSE
  
  n_leaf <- sum(is_leaf)
  
  n_node <- sum(!is_leaf)
  node_id <- ifelse(is_leaf, cumsum(is_leaf), cumsum(!is_leaf) + n_leaf)
  
  # Store the node ids on the graph
  g <- igraph::set_vertex_attr(g, "node_id", index = traverse$order,
                               value = node_id)
  
  # Extract the edge and vertex data
  vertex_data <-  igraph::as_data_frame(g, "vertices") %>%
    mutate(name = row_number())
  edge_data <- igraph::as_data_frame(g, "edges")
  edge_data$geom <- NULL
  
  # Substitute the node id numbers into the edge list
  edge_data <- unlist(edge_data)
  edge_data <- vertex_data$node_id[match(edge_data, vertex_data$name)]
  edge_data <- matrix(edge_data, ncol = 2)
  
  # lookup the tip and node labels
  tip_labels <- 1:n_leaf
  tip_labels <- vertex_data$name[match(tip_labels, vertex_data$node_id)]
  node_labels <- (n_leaf + 1):(n_node + n_leaf)
  node_labels <- vertex_data$name[match(node_labels, vertex_data$node_id)]
  
  
  # Build the phylogeny
  phy <- structure(list(edge = edge_data,
                        edge.length = rep(1, nrow(edge_data)),
                        tip.labels = tip_labels,
                        node.labels = node_labels,
                        Nnode = n_node),
                   class = "phylo")
  
  stra <- as.data.frame(strahlerNumber(phy)) %>%
    rename(strahler_order = `strahlerNumber(phy)`) %>%
    mutate(node_id = row_number()) %>%
    left_join(vertex_data, by = "node_id") %>%
    rename(to = name)
  
  return(stra)
}

ln <- st_read("streams.gpkg") %>%
  st_cast("LINESTRING")

net <- as_sfnetwork(ln)

g <- net %>%
  transposeGraph()

stra <- igraphStrahler(g)

edges = st_as_sf(net, "edges") %>%
  left_join(stra, by = c("from" = "to"))

cols = brewer.pal(3, "Blues")
pal = colorRampPalette(cols)

plot(st_geometry(edges))
plot(edges["strahler_order"], 
     lwd = 3, , 
     add = TRUE,
     col = pal(3)[edges$strahler_order])
legend(x = "topright",
       legend = c("1","2","3"),          
       lwd = 3,
       col = cols,
       title = "Strahler order")