Using sfnetworks for stream network analysis #164

sraul1 · 2021-08-23T21:27:12Z

sraul1
Aug 23, 2021

I posted this on stackoverflow but thought I would also post here as this seems more active of a space.

I'll leave the more specific question I asked on the stackoverflow post if anyone feels like looking at it, but essentially I was just curious of whether anyone had any vignettes/resources for using sfnetworks/tidygraph for hydrologic and stream network analyses purposes. I had not come across much on the subject. Most developers/users of graph analysis tools are more geared toward road network analyses, and as a result most intros on the subject don't deal with river networks. This makes complete sense given the difference in terms of number of users, but there are key differences between these types of applications (namely all stream networks are rooted trees, with only one direction of movement). Given the tidyverse framework that many people in the earth sciences/ecology/GIS are familiar with, it seems like a nature extension to make use of sfnetworks.

These types of network analysis tools would be pretty heavily utilized by people doing hydrology/ecology/general physical science research (which granted, is not a huge population). Especially those of us trying to keep our GIS/analysis/writing in a single space. Most tools available at the moment for hydrologic network analysis are either wrapped up with an ESRI license (All the subcomponents of the Network Analysis toolboxes; ArcHydro), or are 20+ years old written in some variation of FORTRAN and C (SAGA GIS, GRASS GIS, TauDEM).

I'll add that I am a front end user, the folks on the development side are much appreciate by users like myself who less know how to develop packages, and more just know how to Google well enough to get along.

agila5 · 2021-08-25T09:21:46Z

agila5
Aug 25, 2021
Collaborator

Hi @sraul1. First of all, thank you very much for your question. I've never worked with river networks, and that seems a super interesting problem. I tried reading your question on SO and (for me at least) it isn't easy to connect all the pieces. Would you mind simplifying the question on SO starting from a network with just 3 edges and (re)explain the procedure in that case? Something like this:

library(sf)
library(tidygraph)
library(sfnetworks)
library(ggplot2)

n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(2, 1)))
n03 = st_sfc(st_point(c(1, 2)))
n04 = st_sfc(st_point(c(3, 2)))

from = c(1, 2, 2)
to = c(2, 3, 4)

nodes = st_as_sf(c(n01, n02, n03, n04))
edges = data.frame(from = from, to = to)

sfn <- sfnetwork(nodes, edges) %>%
  convert(to_spatial_explicit, .clean = TRUE) %E>%
  mutate(edgeID = c("a","b","c")) %>%
  mutate(dL = c(300, 100, 100))
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

ggplot()+
  geom_sf(
    data = sfn %>% st_as_sf("edges"),
    aes(color = factor(edgeID), linetype = as.character(dL)),
    size = 1.5
  ) +
  geom_sf(data = sfn %>% st_as_sf("nodes"))+
  scale_color_brewer(palette = "Paired")+
  labs(x = "Longitude", y = "Latitude", color = "EdgeID", linetype =  "dL")+
  theme_bw()

^{Created on 2021-08-25 by the reprex package (v2.0.0)}

Sorry for the extra work, but I think that if we solve the simplest case, then it should be easier to generalise the solution to more realistic scenarios. Please also check if the smaller example does make sense.

2 replies

sraul1 Aug 26, 2021
Author

Happy to slim down and hopefully clarify some of the question. I'll try to better described the problem with the 3 edge network procedure, and then again on a slightly more complicated network with 5 edges. I want to add one layer of complexity to example you provided, namely just adding one additional tributary to edge b so that it comprises a true "sub network". Given my starting network is 10,000+ tributaries, the simplification is well justified.

I am working on a solving an equation as you move upstream on a stream, starting at the river outlet. The 2 key parts of the equation that I am trying to use sfnetworks to solve are hx and h0

In gist, h0 needs to be updated as you move upstream to new edges to reflect the maximum hx of the immediate downstream edge.

h0 updates reflect the hx calculated on the most upstream node of the edge immediately downstream. The exception being on the first edge, were h0 is user defined.

The equation and background:

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

N, dX, K, and dL 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 dL 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.
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.
dL 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.
(roughly follows order of solution in graph)

hx is calculated for node 1 of edge a.
hx is calcualted for node 2 of edge a
2.5. Update h0 to reflect the hx calculated on the node 2 of edge a (will be used for edge b and edge c)
hx is calculated for node 1 on edge b
hx is calculated for node 2 on edge b
hx is calculated for node 1 on edge c
hx is calculated for node 2 one edge c

library(sf)
library(tidygraph)
library(sfnetworks)
library(ggplot2)
library(tidyverse)

n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(2, 1)))
n03 = st_sfc(st_point(c(1, 2)))
n04 = st_sfc(st_point(c(3, 2)))

from = c(1, 2, 2)
to = c(2, 3, 4)

nodes = st_as_sf(c(n01, n02, n03, n04))
edges = data.frame(from = from, to = to)

sfn <- sfnetwork(nodes, edges) %>%
  convert(to_spatial_explicit, .clean = TRUE) %E>%
  mutate(edgeID = c("a","b","c")) %>%
  mutate(dL = c(300, 100, 100)) 
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

#easier mapping
df <-
  tibble(
    node = 1:6,
    edgeID = c("a", "a", "c", "c", "b", "b"),
    x = c(0, 2, 2.2, 3, 1.7, 1),
    y = c(0.2, 0.8, 0.9, 1.8, 1.1, 2.1),
    SolutionOrder = c(1, 2, 3, 4, 5, 6),
    N = 0.5,
    K = 1500,
    dX = c(0, 100, 0, 100, 0, 100),
    dL = c(300, 300, 100, 100, 100, 100),
    h0_initial = c(10, 10, NA, NA, NA, NA),
    hx_initial = sqrt(h0_initial ^ 2 + (N * dX / K) * (2 * dL - dX)),
    hx_known = c(10, 10.8, 10.8, 10.95, 10.8, 10.95),
    h0_known = c(10, 10, 10.8, 10.8, 10.8, 10.8),
    h0ValueUsed = c(
      "User defined, 10",
      "User defined, 10",
      "hx calculated at edge a, node 2",
      "hx calculated at edge a, node 2",
      "hx calculated at edge a, node 2",
      "hx calculated at edge a, node 2"
    )
  )

print.data.frame(df)
#>   node edgeID   x   y SolutionOrder   N    K  dX  dL h0_initial hx_initial
#> 1    1      a 0.0 0.2             1 0.5 1500   0 300         10   10.00000
#> 2    2      a 2.0 0.8             2 0.5 1500 100 300         10   10.80123
#> 3    3      c 2.2 0.9             3 0.5 1500   0 100         NA         NA
#> 4    4      c 3.0 1.8             4 0.5 1500 100 100         NA         NA
#> 5    5      b 1.7 1.1             5 0.5 1500   0 100         NA         NA
#> 6    6      b 1.0 2.1             6 0.5 1500 100 100         NA         NA
#>   hx_known h0_known                     h0ValueUsed
#> 1    10.00     10.0                User defined, 10
#> 2    10.80     10.0                User defined, 10
#> 3    10.80     10.8 hx calculated at edge a, node 2
#> 4    10.95     10.8 hx calculated at edge a, node 2
#> 5    10.80     10.8 hx calculated at edge a, node 2
#> 6    10.95     10.8 hx calculated at edge a, node 2
#labels
df2 <- df %>%
  mutate(
    x = c(0, 1.8, 2.3, 3, 1.4, 1.3),
    y = c(0.4, 0.7, 0.8, 1.7, 1.1, 2.1),
    hx = c(10, 10.8, 10.8, 10.95, 10.8, 10.95),
    h0 = c(10, 10, 10.8, 10.8, 10.8, 10.8),
    y1 = y - 0.1,
    y2 = c(3, 7, 13, 21, 10, 17),
    x2 = c(0, 21, 21, 42, 15, 57),
  ) %>% 
  select(node, edgeID, x, y, y1, y2,x2, dX, dL,SolutionOrder, h0, hx)

ggplot()+
  geom_sf(data = sfn %>% st_as_sf("edges"), aes(color = factor(edgeID), linetype = as.character(dL)), size = 1.5) +
  geom_sf(data = sfn %>% st_as_sf("nodes"))+
  scale_color_brewer(palette = "Paired")+
  labs(x = "Longitude", y = "Latitude", color = "EdgeID", linetype =  "dL")+
  theme_bw()+
  geom_label(data = df , aes(x = x, y = y, label = SolutionOrder, color = factor(edgeID)), fontface = "bold",show.legend = F)+
  geom_label(aes(x = 1.5, y = 0.2, label = "Numbers indicate order the nodes are solved"))+
  geom_text(data = df2, aes(x = x, y = y, label = paste("h0 =", sep = " ", h0), color = factor(edgeID)),show.legend = F)+
  geom_text(data = df2, aes(x = x, y = y1, label = paste("hx =", sep = " ", hx), color = factor(edgeID)),show.legend = F)

^{Created on 2021-08-26 by the reprex package (v2.0.1)}

The problem in a slightly more complicated network with 5 edges:

Start the calculation at the downstream node of edge a with h0 = 10.
(roughly follows order of solution in graph)

hx is calculated for node 1 of edge a.
hx is calcualted for node 2 of edge a
2.5. Update h0 to reflect the hx calculated on the node 2 of edge a (will be used for edge b and edge c)
hx is calculated for node 1 on edge b
hx is calculated for node 2 on edge b
hx is calculated for node 1 on edge c
hx is calculated for node 2 one edge c
6.5. Update h0 to reflect the hx calculated on node 2 of edge c (will be used for edge d and edge e)
hx is calculated for node 1 of edge d
hx is calculated for node 2 of edge d
hx is calculated for node 1 of edge e
hx is calculated for node 2 of edge e

library(sf)
library(tidygraph)
library(sfnetworks)
library(ggplot2)
library(tidyverse)

n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(2, 1)))
n03 = st_sfc(st_point(c(3, 2)))
n04 = st_sfc(st_point(c(1.5, 1.5)))
n05 = st_sfc(st_point(c(2, 2)))
n06 = st_sfc(st_point(c(1, 2)))

from = c(1, 2, 2, 4, 4)
to = c(2, 3, 4, 5, 6)

nodes = st_as_sf(c(n01, n02, n03, n04, n05, n06))
edges = data.frame(from = from, to = to)

sfn <- sfnetwork(nodes, edges) %>%
  convert(to_spatial_explicit, .clean = TRUE) %E>%
  mutate(edgeID = c("a","b","c", "d", "e")) %>%
  mutate(dL = c(500, 100, 300, 100, 100))
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

#easier mapping
df <- tibble(
  x = c(0, 2, 2.2, 2, 3, 1.5, 1.5, 1.3, 2.2, 0.8),
  y = c(0.1, 0.8, 1, 1.2, 1.8, 1.3, 1.7, 1.5, 2, 2),
  node = c(1, 1, 2, 2, 3, 4, 4, 4, 5, 6),
  edgeID = c("a", "a", "b", "c", "b", "c", "d", "e", "d", "e"),
  SolutionOrder = c(1, 2, 3, 5, 4, 6, 7, 9, 8, 10),
  N = 0.5,
  K = 1500,
  dX = c(0, 100, 0, 0, 100, 100, 0, 0, 100, 100),
  dL = c(500, 500, 100, 300, 100, 300, 100, 100, 100, 100),
  h0_initial = c(10, 10, NA, NA, NA, NA, NA, NA, NA, NA),
  hx_initial = sqrt(h0_initial ^ 2 + (N * dX / K) * (2 * dL - dX)),
  hx_known = c(10, 11.4, 11.4, 11.4, 11.54, 12.1, 12.1, 12.1, 12.24, 12.24),
  h0_known = c(10, 10, 11.4, 11.4, 11.4, 11.4, 12.1, 12.1, 12.1, 12.1),
  h0ValueUsed = c(
    "User defined, 10",
    "User defined, 10",
    "hx calculated at edge a, node 2",
    "hx calculated at edge a, node 2",
    "hx calculated at edge a, node 2",
    "hx calculated at edge a, node 2",
    "hx calculated at edge c, node 4",
    "hx calculated at edge c, node 4",
    "hx calculated at edge c, node 4",
    "hx calculated at edge c, node 4"
  )
)

print.data.frame(df)
#>      x   y node edgeID SolutionOrder   N    K  dX  dL h0_initial hx_initial
#> 1  0.0 0.1    1      a             1 0.5 1500   0 500         10   10.00000
#> 2  2.0 0.8    1      a             2 0.5 1500 100 500         10   11.40175
#> 3  2.2 1.0    2      b             3 0.5 1500   0 100         NA         NA
#> 4  2.0 1.2    2      c             5 0.5 1500   0 300         NA         NA
#> 5  3.0 1.8    3      b             4 0.5 1500 100 100         NA         NA
#> 6  1.5 1.3    4      c             6 0.5 1500 100 300         NA         NA
#> 7  1.5 1.7    4      d             7 0.5 1500   0 100         NA         NA
#> 8  1.3 1.5    4      e             9 0.5 1500   0 100         NA         NA
#> 9  2.2 2.0    5      d             8 0.5 1500 100 100         NA         NA
#> 10 0.8 2.0    6      e            10 0.5 1500 100 100         NA         NA
#>    hx_known h0_known                     h0ValueUsed
#> 1     10.00     10.0                User defined, 10
#> 2     11.40     10.0                User defined, 10
#> 3     11.40     11.4 hx calculated at edge a, node 2
#> 4     11.40     11.4 hx calculated at edge a, node 2
#> 5     11.54     11.4 hx calculated at edge a, node 2
#> 6     12.10     11.4 hx calculated at edge a, node 2
#> 7     12.10     12.1 hx calculated at edge c, node 4
#> 8     12.10     12.1 hx calculated at edge c, node 4
#> 9     12.24     12.1 hx calculated at edge c, node 4
#> 10    12.24     12.1 hx calculated at edge c, node 4

#labels
df2 <- df %>% 
  mutate(
    x = c(0, 2,2.5,2,2.9,1.5,1.5,1,2.5,0.8),
    y = c(0.3,0.7,1,1.45,1.7,1.2,1.9,1.5,2,1.9),
    hx = c(10, 11.4, 11.4, 11.4, 11.54, 12.1, 12.1, 12.1, 12.24, 12.24),
    h0 = c(10, 10, 11.4, 11.4, 11.4, 11.4, 12.1, 12.1, 12.1, 12.1),
    y1 = y - 0.1,
    y2 = c(0.3,0.8,1.2,1.3,1.7,1.2,1.7,1.5,1.7,1.7),
    x2 = c(0, 2,2.3,2,3,1.5,1.7,1.3,2.3,0.8),
  ) %>% 
  select(node, edgeID, x, y, y1, y2,x2, dX, dL,SolutionOrder, h0, hx)
  
ggplot()+
  geom_sf(data = sfn %>% st_as_sf("edges"),aes(color = factor(edgeID), linetype = as.character(dL)),size = 1.5) +
  geom_sf(data = sfn %>% st_as_sf("nodes"))+
  scale_color_brewer(palette = "Paired")+
  scale_linetype_manual(values = c("dotted", "twodash", "solid"))+
  labs(x = "Longitude", y = "Latitude", color = "EdgeID", linetype =  "dL")+
  theme_bw()+
  geom_label(data = df , aes(x = x, y = y, label = SolutionOrder, color = factor(edgeID)), fontface = "bold",show.legend = F)+
  geom_label(aes(x = 1.5, y = 0.2, label = "Numbers indicate order the nodes are solved"))+
  geom_text(data = df2, aes(x = x, y = y, label = paste("h0 =", sep = " ", h0), color = factor(edgeID)),show.legend = F)+
  geom_text(data = df2, aes(x = x, y = y1, label = paste("hx =", sep = " ", hx), color = factor(edgeID)),show.legend = F)

^{Created on 2021-08-26 by the reprex package (v2.0.1)}

In trying to use some language from graph theory, I was hoping to use sfnetworks to retain the parent-child relationship between nodes and edges as a way to iterate through the entire stream network solving this equation. After the initial downstream calculation at edge a, node 1, every child needs to know the maximum hx of it's parent before being able to calculate hx for itself.

Happy to clarify more, I realize I may not be presenting the question clearly.

agila5 Aug 26, 2021
Collaborator

Dear @sraul1, thank you very much for the new examples. They are simply amazing. I'll try to answer your question sometime during the next week and (hopefully) add some more comments here as soon as possible.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Using sfnetworks for stream network analysis #164

{{title}}

Replies: 1 comment 2 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Using sfnetworks for stream network analysis #164

sraul1 Aug 23, 2021

Replies: 1 comment · 2 replies

agila5 Aug 25, 2021 Collaborator

sraul1 Aug 26, 2021 Author

agila5 Aug 26, 2021 Collaborator

sraul1
Aug 23, 2021

Replies: 1 comment 2 replies

agila5
Aug 25, 2021
Collaborator

sraul1 Aug 26, 2021
Author

agila5 Aug 26, 2021
Collaborator