1

I have a dataset that contains network distances along a river network. This dataset was created using Topotoolbox in Matlab, with the node and segment information exported. I've been working on the analysis in R mainly due to familiarity, but am open to just about any other software (GIS or otherwise).

I have 5 river segments, segments 1-5. Each segment has 5 nodes on it (index). The connection between nodes are dictated by the ix (giver) and ixc (receiver) value. The distance along each segment moving upstream is stored in seg.dist (0 at the outlet, 40 at the top), and the distance from the network outlet, or root, stored in net.dist. The total distance upstream of each segment is stored in up.dist. y, z just contains values for mapping. n and k are constants.

library(tidyverse)
segment <- rep(1:5, each = 5)
index <- rep(1:5, times = 5)
ix <- 1:25
ixc <- c(0,1,2,3,4,5,6,7,8,9,5,11,12,13,14,10,16,17,18,19,10,21,22,23,24)
seg.dist <- rep(c(0,10,20,30,40), times = 5)
net.dist <- c(0,10,20,30,40,40,50,60,70,80,40,50,60,70,80,80,90,100,110,120,80,90,100,110,120)
up.dist <- c(400,400,400,400,400,200,200,200,200,200,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
y <-c(0,0,0,0,0,0,1,1,1,1,0,-1.2,-1.4,-1.6,-1.8,1,1.2,1.4,1.6,1.8,1,0.8,0.6,0.4,0.2)
z <- c(1,2,3,4,5,5,6,7,8,9,5,7,9,11,13,9,10,11,12,13,9,11,13,15,17)
df <- data.frame(segment, index, ix, ixc, seg.dist, net.dist,up.dist,y, z) %>% 
  mutate(segment = factor(segment)) %>% 
  mutate(n = 0.5, k = 1500)
> df
   segment index ix ixc seg.dist net.dist up.dist    y  z   n    k
1        1     1  1   0        0        0     400  0.0  1 0.5 1500
2        1     2  2   1       10       10     400  0.0  2 0.5 1500
3        1     3  3   2       20       20     400  0.0  3 0.5 1500
4        1     4  4   3       30       30     400  0.0  4 0.5 1500
5        1     5  5   4       40       40     400  0.0  5 0.5 1500
6        2     1  6   5        0       40     200  0.0  5 0.5 1500
7        2     2  7   6       10       50     200  1.0  6 0.5 1500
8        2     3  8   7       20       60     200  1.0  7 0.5 1500
9        2     4  9   8       30       70     200  1.0  8 0.5 1500
10       2     5 10   9       40       80     200  1.0  9 0.5 1500

I want to iteratively make a calculation of a variable called h, starting from the bottom-most segment, and index (segment 1, index 1), moving upward through the index and the segment value. h is calculated using the variables defined above, but with an additional variable called h0.

h = sqrt(h0^2 + ((n*seg.dist)/k)*(2*up.dist - seg.dist))

h0 serves as a sort of initialization variable, that I want to iteratively update as h is calculated. At this lowest segment, I want to initialize h0 to a value of 10.

#Segment 1
h0 = 10
h1.1 = sqrt(h0^2 + (0.5*0/1500)*(2*400 - 0))
h1.2 = sqrt(h0^2 + (0.5*10/1500)*(2*400 - 10))
h1.3 = sqrt(h0^2 + (0.5*20/1500)*(2*400 - 20))
h1.4 = sqrt(h0^2 + (0.5*30/1500)*(2*400 - 30))
h1.5 = sqrt(h0^2 + (0.5*40/1500)*(2*400 - 40))
> h1.1
[1] 10
> h1.2
[1] 10.13081
> h1.3
[1] 10.25671
> h1.4
[1] 10.37786
> h1.5
[1] 10.49444

When I get to the segment 2 and segment 3, h0 changes. I want to update it so that:

h0 = max(h(parent segment)).

So, for segments 2 and 3, h0 = h1.5.

h0 = h1.5 
#Segment 2
h2.1 = sqrt(h0^2 + (0.5*0/1500)*(2*200 - 0))
h2.2 = sqrt(h0^2 + (0.5*10/1500)*(2*200 - 10))
h2.3 = sqrt(h0^2 + (0.5*20/1500)*(2*200 - 20))
h2.4 = sqrt(h0^2 + (0.5*30/1500)*(2*200 - 30))
h2.5 = sqrt(h0^2 + (0.5*40/1500)*(2*200 - 40))
#Segment 3
h3.1 = sqrt(h0^2 + (0.5*0/1500)*(2*0 - 0))
h3.2 = sqrt(h0^2 + (0.5*10/1500)*(2*0 - 10))
h3.3 = sqrt(h0^2 + (0.5*20/1500)*(2*0 - 20))
h3.4 = sqrt(h0^2 + (0.5*30/1500)*(2*0 - 30))
h3.5 = sqrt(h0^2 + (0.5*40/1500)*(2*0 - 40))

Again, when we get to segments 4 and 5, h0 changes. It will be equal to the maximum h calculation from segment 2.

> h2.1
[1] 10.49444
> h2.2
[1] 10.5562
> h2.3
[1] 10.61446
> h2.4
[1] 10.66927
> h2.5
[1] 10.7207
h0 = h2.5 
#Segment 4
h4.1 = sqrt(h0^2 + (0.5*0/1500)*(2*0 - 0))
h4.2 = sqrt(h0^2 + (0.5*10/1500)*(2*0 - 10))
h4.3 = sqrt(h0^2 + (0.5*20/1500)*(2*0 - 20))
h4.4 = sqrt(h0^2 + (0.5*30/1500)*(2*0 - 30))
h4.5 = sqrt(h0^2 + (0.5*40/1500)*(2*0 - 40))
#Segment 5
h5.1 = sqrt(h0^2 + (0.5*0/1500)*(2*0 - 0))
h5.2 = sqrt(h0^2 + (0.5*10/1500)*(2*0 - 10))
h5.3 = sqrt(h0^2 + (0.5*20/1500)*(2*0 - 20))
h5.4 = sqrt(h0^2 + (0.5*30/1500)*(2*0 - 30))
h5.5 = sqrt(h0^2 + (0.5*40/1500)*(2*0 - 40))

The final data.frame looks something like this:

h.df <- data.frame(h1.1,h1.2,h1.3,h1.4,h1.5,h2.1,h2.2,h2.3,h2.4,h2.5,h3.1,h3.2,h3.3,h3.4,h3.5,h4.1,h4.2,h4.3,h4.4,h4.5,h5.1,h5.2,h5.3,h5.4,h5.5) %>% 
  pivot_longer(everything(),  names_to = NULL, values_to = "h")

df.final <- bind_cols(df, h.df) %>% 
  mutate(h0 = c(10,10,10,10,10,10.49,10.49,10.49,10.49,10.49,10.49,10.49,10.49,10.49,10.49, 10.72, 10.72, 10.72, 10.72, 10.72, 10.72, 10.72, 10.72, 10.72, 10.72))

> df.final
   segment index ix ixc seg.dist net.dist up.dist    y  z   n    k        h    h0
1        1     1  1   0        0        0     400  0.0  1 0.5 1500 10.00000 10.00
2        1     2  2   1       10       10     400  0.0  2 0.5 1500 10.13081 10.00
3        1     3  3   2       20       20     400  0.0  3 0.5 1500 10.25671 10.00
4        1     4  4   3       30       30     400  0.0  4 0.5 1500 10.37786 10.00
5        1     5  5   4       40       40     400  0.0  5 0.5 1500 10.49444 10.00
6        2     1  6   5        0       40     200  0.0  5 0.5 1500 10.49444 10.49
7        2     2  7   6       10       50     200  1.0  6 0.5 1500 10.55620 10.49
8        2     3  8   7       20       60     200  1.0  7 0.5 1500 10.61446 10.49
9        2     4  9   8       30       70     200  1.0  8 0.5 1500 10.66927 10.49
10       2     5 10   9       40       80     200  1.0  9 0.5 1500 10.72070 10.49
11       3     1 11   5        0       40       0  0.0  5 0.5 1500 10.49444 10.49
12       3     2 12  11       10       50       0 -1.2  7 0.5 1500 10.49285 10.49
13       3     3 13  12       20       60       0 -1.4  9 0.5 1500 10.48809 10.49
14       3     4 14  13       30       70       0 -1.6 11 0.5 1500 10.48014 10.49
15       3     5 15  14       40       80       0 -1.8 13 0.5 1500 10.46900 10.49
16       4     1 16  10        0       80       0  1.0  9 0.5 1500 10.72070 10.72
17       4     2 17  16       10       90       0  1.2 10 0.5 1500 10.71914 10.72
18       4     3 18  17       20      100       0  1.4 11 0.5 1500 10.71448 10.72
19       4     4 19  18       30      110       0  1.6 12 0.5 1500 10.70670 10.72
20       4     5 20  19       40      120       0  1.8 13 0.5 1500 10.69579 10.72
21       5     1 21  10        0       80       0  1.0  9 0.5 1500 10.72070 10.72
22       5     2 22  21       10       90       0  0.8 11 0.5 1500 10.71914 10.72
23       5     3 23  22       20      100       0  0.6 13 0.5 1500 10.71448 10.72
24       5     4 24  23       30      110       0  0.4 15 0.5 1500 10.70670 10.72
25       5     5 25  24       40      120       0  0.2 17 0.5 1500 10.69579 10.72

I have made some runs at a solution for this using graph theory approaches but have had limited success. The problem seems recursive in nature, which I have limited experience with. It seems as though a structure of parent/child relationships is going to be required in solving this, similar to a tree data structure or list of list. I have presented this problem as simply as I could, but the reality is the solution will be applied to 10k+ segment, with hundreds of nodes (index) within each segment. This kind of problem seems like it would frequently come up in many network/graph theory work.

UPDATE

@gavg712 your answer works great for the example I presented, I have been attempting to alter it to more realistic networks. The linear referencing of parent/child isn't the best so I have working to update your solution. This is based on @Spacedman comment of using a parentSegment column.

Say I add a column named parent to the original df, this column reflects the segment ID of the parent. If a segment has no parent then it has the value of 0 (the starting point for the sequence/calculation).

df <- df %>% 
  mutate(parent = c(0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2))

> df
   segment index ix ixc seg.dist net.dist up.dist    y  z   n    k parent
1        1     1  1   0        0        0     400  0.0  1 0.5 1500      0
2        1     2  2   1       10       10     400  0.0  2 0.5 1500      0
3        1     3  3   2       20       20     400  0.0  3 0.5 1500      0
4        1     4  4   3       30       30     400  0.0  4 0.5 1500      0
5        1     5  5   4       40       40     400  0.0  5 0.5 1500      0
6        2     1  6   5        0       40     200  0.0  5 0.5 1500      1
7        2     2  7   6       10       50     200  1.0  6 0.5 1500      1
8        2     3  8   7       20       60     200  1.0  7 0.5 1500      1
9        2     4  9   8       30       70     200  1.0  8 0.5 1500      1
10       2     5 10   9       40       80     200  1.0  9 0.5 1500      1
11       3     1 11   5        0       40       0  0.0  5 0.5 1500      1
12       3     2 12  11       10       50       0 -1.2  7 0.5 1500      1
13       3     3 13  12       20       60       0 -1.4  9 0.5 1500      1
14       3     4 14  13       30       70       0 -1.6 11 0.5 1500      1
15       3     5 15  14       40       80       0 -1.8 13 0.5 1500      1
16       4     1 16  10        0       80       0  1.0  9 0.5 1500      2
17       4     2 17  16       10       90       0  1.2 10 0.5 1500      2
18       4     3 18  17       20      100       0  1.4 11 0.5 1500      2
19       4     4 19  18       30      110       0  1.6 12 0.5 1500      2
20       4     5 20  19       40      120       0  1.8 13 0.5 1500      2
21       5     1 21  10        0       80       0  1.0  9 0.5 1500      2
22       5     2 22  21       10       90       0  0.8 11 0.5 1500      2
23       5     3 23  22       20      100       0  0.6 13 0.5 1500      2
24       5     4 24  23       30      110       0  0.4 15 0.5 1500      2
25       5     5 25  24       40      120       0  0.2 17 0.5 1500      2

How would you alter the function so that instead of updating h0 based on the previous segment maximum (based on a linear indexing reference system between segments), it instead uses that parent column as the index, and updates h0 as the maximum from the parent segment. Something along the lines of:

update <- function(seg, h0, n, seg.dist, k, up.dist, parent, output = c("h", "h0")){
  ctrl <- c(0, diff(seg))
  h <- list(h = numeric(), h0 = numeric())
  j = 1
  for(i in seq_along(seg)) {
    if(ctrl[i] > 0) {
      h0 = max(h[[1]][parent])
      j = i
    }
    h[[1]][i] <- sqrt(h0^2 + ((n[i]*seg.dist[i])/k[i])*(2*up.dist[i] - seg.dist[i]))
    h[[2]][i] <- h0
  }
  if(length(output) == 1) return(h[[output]]) else return(h[output]) 
}

Hopefully that is clear, this type of solution does seem to be teetering the edge of needing a more involved recursive solution in order to be applied to my full dataset.

2
  • 1
    h for segment block i is dependent only on one segment block j < i, so you can start with an NA column of h and fill each segment block in increasing i. A column of parentSegment might be useful. I'm not sure if the real data is also going to be ordered like this, but a pre-processing step based on distance from the head node would make it so.
    – Spacedman
    Feb 12 at 21:25
  • Your parentSegment suggestion is exactly the direction I am trying to move towards. @gavg712 solution works great for the simple example presented here. The linear indexing of the reference system is a bit of an issue on the actual dataset when you have more than one segment come to a confluence. Segments 2 and 3 may come together at segment 1, but the solution at the moment would try to use the max h from segment 2 as the h0 for segment 3, instead of the max h from segment 1. The segments are labelled from 1 - n, moving upstream.
    – Scott R
    Mar 5 at 13:22
3

Hope is what you are looking for:


# function to calculate both
h <- function(seg, h0 = 10, n, seg.dist, k, up.dist, output = c("h", "h0")){
  ctrl <- c(0, diff(seg))
  h <- list(h = numeric(), h0 = numeric())
  j = 1
  for(i in seq_along(seg)) {
    if(ctrl[i] > 0) {
      h0 = max(h[[1]][j:(i-1)])
      j = i
    }
    h[[1]][i] <- sqrt(h0^2 + ((n[i]*seg.dist[i])/k[i])*(2*up.dist[i] - seg.dist[i]))
    h[[2]][i] <- h0
  }
  if(length(output) == 1) return(h[[output]]) else return(h[output]) 
}

# testing function with name of variable
df %>%
  mutate(h = h(segment, h0 = 10, n, seg.dist, k, up.dist, "h"),
         h0 = h(segment, h0 = 10, n, seg.dist, k, up.dist, "h0"))
#>    segment index ix ixc seg.dist net.dist up.dist    y  z   n    k        h
#> 1        1     1  1   0        0        0     400  0.0  1 0.5 1500 10.00000
#> 2        1     2  2   1       10       10     400  0.0  2 0.5 1500 10.13081
#> 3        1     3  3   2       20       20     400  0.0  3 0.5 1500 10.25671
#> 4        1     4  4   3       30       30     400  0.0  4 0.5 1500 10.37786
#> 5        1     5  5   4       40       40     400  0.0  5 0.5 1500 10.49444
#> 6        2     1  6   5        0       40     200  0.0  5 0.5 1500 10.49444
#> 7        2     2  7   6       10       50     200  1.0  6 0.5 1500 10.55620
#> 8        2     3  8   7       20       60     200  1.0  7 0.5 1500 10.61446
#> 9        2     4  9   8       30       70     200  1.0  8 0.5 1500 10.66927
#> 10       2     5 10   9       40       80     200  1.0  9 0.5 1500 10.72070
#> 11       3     1 11   5        0       40       0  0.0  5 0.5 1500 10.72070
#> 12       3     2 12  11       10       50       0 -1.2  7 0.5 1500 10.71914
#> 13       3     3 13  12       20       60       0 -1.4  9 0.5 1500 10.71448
#> 14       3     4 14  13       30       70       0 -1.6 11 0.5 1500 10.70670
#> 15       3     5 15  14       40       80       0 -1.8 13 0.5 1500 10.69579
#> 16       4     1 16  10        0       80       0  1.0  9 0.5 1500 10.72070
#> 17       4     2 17  16       10       90       0  1.2 10 0.5 1500 10.71914
#> 18       4     3 18  17       20      100       0  1.4 11 0.5 1500 10.71448
#> 19       4     4 19  18       30      110       0  1.6 12 0.5 1500 10.70670
#> 20       4     5 20  19       40      120       0  1.8 13 0.5 1500 10.69579
#> 21       5     1 21  10        0       80       0  1.0  9 0.5 1500 10.72070
#> 22       5     2 22  21       10       90       0  0.8 11 0.5 1500 10.71914
#> 23       5     3 23  22       20      100       0  0.6 13 0.5 1500 10.71448
#> 24       5     4 24  23       30      110       0  0.4 15 0.5 1500 10.70670
#> 25       5     5 25  24       40      120       0  0.2 17 0.5 1500 10.69579
#>          h0
#> 1  10.00000
#> 2  10.00000
#> 3  10.00000
#> 4  10.00000
#> 5  10.00000
#> 6  10.49444
#> 7  10.49444
#> 8  10.49444
#> 9  10.49444
#> 10 10.49444
#> 11 10.72070
#> 12 10.72070
#> 13 10.72070
#> 14 10.72070
#> 15 10.72070
#> 16 10.72070
#> 17 10.72070
#> 18 10.72070
#> 19 10.72070
#> 20 10.72070
#> 21 10.72070
#> 22 10.72070
#> 23 10.72070
#> 24 10.72070
#> 25 10.72070

Created on 2021-02-15 by the reprex package (v1.0.0)

Best

1
  • This was very helpful, I updated the post based on some of the progress made with this potential solution.
    – Scott R
    Mar 3 at 22:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.