Последовательность расчета речной сети с использованием sfnetworks и r

Question

У меня есть корневое дерево. Он получен из шейп-файла линии потоковой сети. Для каждого притока в потоковой сети, очевидно, есть начальный и конечный узлы, которые определяют, но есть также узлы, соединяющие начальный и конечный узлы на линии. Для простоты я возьму пример сети, разработанной в этом вопросе .

      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() 

Сетевой график(почему-то у меня нет 10 репутации в обычном stackoverflow, извиняюсь за картинку)

Я работаю над решением уравнения, когда вы движетесь вверх по течению, начиная с устья реки. Я сохраню вам как можно больше деталей этого уравнения, но по сути:

Я решаю для каждого узла x на основе атрибутов узла.

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

• `` и - все фиксированные параметры, известные до начала расчета.

• , фиксированы для всех узлов, они никогда не меняются.

• является специфическим атрибутом узла.

• и являются атрибутами, специфичными для ребра (они одинаковы для всех узлов на ребре)

• определяется пользователем в начале расчета, но затем изменяется по мере прохождения сети в восходящем направлении, изменяясь в конечном узле, когда сходятся новые ребра.

• изначально установлено на 10, N = 0,5, K = 1500.

• Для простоты изложения мы скажем, что каждое ребро имеет длину 100 м.

• просто отражает расстояние на каждом краю. Нисходящий узел на каждом краю, . Восходящий узел на каждом краю, .

• отражает общее расстояние канала до канала плюс общую длину этого канала. .

Я хочу начать вычисление на нижнем узле ребра a с .

1. вычисляется в этой начальной точке, а затем в каждом узле вдоль кромки (в этом примере показаны только начальный и конечный узел), пока вы не дойдете до соединения, где кромка g и кромка b сливаются, образуя кромку a .

2. В этой конечной точке ребра a :

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

3. На этом стыке ребра g и ребра b я хочу обновить, чтобы отразить вычисленное в конечной точке ребра a , поэтому 13.03. Это значение будет использоваться как для кромки g, так и для кромки b .

4. Теперь я хочу сделать этот расчет для каждой подсети. Я начну с подсети, которая начинается на ребре b . был изменен на 13.03, мы вычисляем для узлов на ребре b , затем переходим к стыку ребра c и ребра d .

5. В этой конечной точке ребра b :

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

6. обновляется как для ребра c, так и для ребра d, чтобы он отражал конечную точку ребра b (14.13)
7. вычисляется для всех узлов на ребре d и снова на ребре c .
8. В конечной точке ребра c обновляется, чтобы отразить точку соединения ребра e и ребра f.

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

9. окончательно вычисляется для всех узлов на ребре e и f , используя значение, вычисленное в конечной точке ребра c для. На этом обход этой подсети завершен.

10. Мы возвращаемся к ребру g, чтобы вычислить подсеть, которая там начинается.

11. Как мы знаем из предыдущих расчетов, в конечной точке ребра а равно 13,03. Теперь мы хотим снова отразить это значение для начала обхода подсети, начиная с края g .

12. Когда мы подходим к стыку узлов на ребрах h и i , мы снова обновляем, чтобы отразить вычисленное в конечной точке ребра g .

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

13. Это значение, используемое для ребер h и i . Мы вычисляем по ребру i, а затем по ребру h . обновляется в последний момент, когда кромка j и кромка k сливаются, чтобы отразить вычисленное в конечной точке кромки h .

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

После этого окончательного обновления вычисляется для всех узлов на ребре j и ребре k .

По сути, необходимо обновлять по мере продвижения вверх по потоку к новым краям, чтобы отразить максимум непосредственного нижнего края. Существует определенный порядок решения, который должен быть достигнут в том смысле, что перед решением любого ребра (за пределами первого ребра) вам необходимо знать для этого сегмента, то есть сначала нужно решить край ниже, чтобы определить использовать.

Я изо всех сил пытался придумать решение цикла / итерации / рекурсии для этого вычисления, которое на первый взгляд кажется, что это должно быть довольно легко решено. Решение по единому дренажному потоку очень легко сделать вручную, но когда оно масштабируется до всей сети, возникают эти проблемы с последовательностью, которые необходимо решить. Когда у вас много тысяч ребер в сети, выполнение этих вычислений вручную невозможно.

Этот тип последовательности расчетов довольно распространен во многих основных гидрологических анализах (расчет порядка потока, величины потока, расчет площади, вносящей вклад вверх по течению, и т. Д.). и кажутся хорошо подходящими для выполнения этого типа вычислений, но есть несколько примеров приложений для гидрологических сетей (большинство инструментов для анализа графиков ориентированы на дорожные сети, что, учитывая разницу в количестве пользователей, имеет смысл).

Я постарался сделать это как можно больше репексом, если потребуется больше, я с радостью предоставлю это, это было сложно сделать, учитывая более концептуальный характер вопроса.