Strahler stream order using igraph or sfnetwork in R

[jahjones]

I cannot fathom how to derive strahler order in R for large networks. Here's an example in postgres and neo4j. I have a solution that works for small networks.

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 a solution that relies on phytools::strahlerNumber which I had to modify. But it doesn't work for larger networks.

e.g. a network with 28 nodes runs in seconds and a network with 1,467 nodes causes R to abort.

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")

[agila5]

Hi @jahjones and thanks for your Q. I'll try to check your problem as soon as possible (hopefully next week). Meanwhile, I think you may want to read the solution here for inspiration: https://stackoverflow.com/questions/65958776/group-id-for-all-edges-between-points-using-sfnetwork-in-r

Moreover, what do you mean by "it doesn't work for larger networks"? Does it take too much time? Are you running out of RAM? Something else?

[jahjones]

e.g. a network with 28 nodes runs in seconds and a network with 1,467 nodes causes R to abort. The offending function in phytools::strahlerNumber() is pw<-reorder(tree,"pruningwise").