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 blocki
is dependent only on one segment blockj < i
, so you can start with anNA
column ofh
and fill each segment block in increasingi
. A column ofparentSegment
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.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 onesegment
come to a confluence.Segments
2 and 3 may come together atsegment
1, but the solution at the moment would try to use the maxh
fromsegment
2 as theh0
forsegment
3, instead of the maxh
fromsegment
1. Thesegments
are labelled from 1 - n, moving upstream.