Sequencing of river network calculation using sfnetworks and r

Question

I’ve got a rooted tree setup in sfnetworks. It is derived from a line shapefile of a stream network. For each tributary in the stream network, there are obviously start and end nodes which sfnetworks determines, but there are also nodes connecting between start and end nodes on the line. I am going to take the example network developed in this question for simplicity sake.

library(sfnetworks)
library(sf)
library(tidygraph)
library(dplyr)
library(tidyverse)
rm(list = ls())

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_1 = sfnetwork(nodes, edges)%>%
  convert(to_spatial_explicit, .clean = TRUE) %>% 
  activate(edges) %>% 
  mutate(edgeID = c("a","c","d","e","f","b","g","h","i","j","k")) %>%
  mutate(dL = c(1100, 300, 100, 100, 100, 500, 500, 300, 100, 100, 100))  


ggplot()+
  geom_sf(data = G_1 %>% activate(edges) %>% as_tibble() %>% st_as_sf(), aes(color = factor(edgeID)), size = 1.5)+
  geom_sf(data = G_1 %>% activate(nodes) %>% as_tibble() %>% st_as_sf())+
  scale_color_brewer(palette = "Paired")+
  labs(x = "Longitude", y = "Latitude", color = "EdgeID")+
  theme_bw() 

I am working on a solving an equation as you move upstream on a stream, starting at the river outlet. I'll save you as much of the details on that equation as I can, but essentially:

I am solving for hx at each node, x, based on node attributes.

hx = sqrt(h0^2 + (N*dX/K)*(2*L – dX))

• N, dX, K, and L are all fixed parameters known before we start the calculation.

• N, K are fixed for all nodes, they never change.

• dX is a node specific attribute.

• h0 and L are edge specific attributes (they are the same for all nodes on an edge)

• h0 is defined by the user at the start of the calculation, but then changes as you traverse the network upstream, changing at the end node when new edges converge.

• h0 is initially set to 10, N = 0.5, K = 1500.

• For sake of simplicity of presenting this, we will say that each edge is 100m long.

• dX just reflects the distance on each edge. The downstream node on each edge, dX = 0. The upstream node on each edge, dX = 100.

• L reflects the total distance of channel upstream of an edge, plus the total length of that edge.

I want to start the calculation at the downstream node of edge a with h0 = 10.

1. hx is calculated at that starting point and then at every node along the edge (in this example only the starting and ending node is shown) until you get to the junction where edge g and edge b merge to form edge a.

2. At that end point of edge a:

hx = sqrt(10^2 + (0.5*100/1500)*(2*1100-100)).
hx = 13.03

3. At this junction of edge g and edge b, I want to update h0 to reflect the hx calculated at the end point of edge a, so 13.03. This hx value will be used as h0 for both edge g and edge b.

4. I now want to make this calculation for each subnetwork. I will start with the subnetwork that starts at edge b. h0 has been changed to 13.03, we calculate hx for the nodes on edge b, then get to the junction of edge c and edge d.

5. At that end point of edge b:

hx = sqrt(13.03^2 + (0.5*100/1500)*(2*500-100))
hx = 14.13

6. h0 is updated for both edge c and edge d so that it reflects hx of edge b’s endpoint (14.13)
7. hx is calculated for all nodes on edge d, and again on edge c.
8. At the endpoint of edge c, h0 is updated to reflect the hx at that junction of edge e and edge f

hx = sqrt(14.13^2+(0.5*100/1500)*(2*300 – 100))
hx = 14.79

9. hx is finally calculated for all the nodes on edge e and edge f, using that hx calculated at the endpoint of edge c for the h0. That completes the traversal of that subnetwork.

10. We return to edge g to calculate hx for the subnetwork that begins there.

11. As we know from previous calculations, hx at the endpoint of edge a is 13.03. We now want h0 to again reflect this value for the beginning of the traversal of the subnetwork starting at edge g.

12. When we get to the node junction at edges h and edge i, again, we update h0 to reflect the hx calculated at the end point of edge g.

hx = sqrt(13.03^2 +(0.5*100/1500)*(2*500-100))
hx = 14.13

13. This is the value used as h0 for edges h and edge i. We calculate hx on edge i and then on edge h. h0 is updated a final time where edge j and edge k merge to reflect the hx calculated at the endpoint of edge h.

hx = sqrt(14.13^2 + (0.5*100/1500)*(2*500-100))
hx = 14.79

After that final update to h0, hx is calculated for all nodes on edge j and edge k.

In gist, h0 needs to be updated as you move upstream to new edges to reflect the maximum hx of the immediate downstream edge. There is a certain order of solution that must be achieved in that, before solving any edge (outside of the first edge), you need to know the ho for that segment, meaning you have to solve the edge below first, to determine the hx to use.

I have struggled to come up with a loop/iteration/recursion solution to this calculation that at the surface seems like it should be quite easily solved. The solution along a single drainage thread is very easy to do manually, but when it is scaled up to a whole network there are these sequencing concerns that must be addressed. When you have many thousand of edges in a network, manually making these calculations would not be an option.

This type of calculation sequencing is quite common in many basic hydrologic analyses (calculating stream order, stream magnitude, calculating upstream contributing area, etc.). sfnetworks and tidygraph seem well suited to make this type of calculation but there are few example applications on hydrologic networks (most tools for graph analysis are geared for road networks which, given the difference in terms of number of users, makes complete sense).

I tried to make this as much of a repex as possible, if more is needed I will happily supply that, it was difficult to do given the more conceptual nature of the question.