# Recursive calculations on a river network

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.

• `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. Feb 12, 2021 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. Mar 5, 2021 at 13:22

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

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