1

我有一个带有空间显式边缘 ( ln_sfnetwork) 的有根树，以及通过添加点层 ( ) 创建的附加边缘pt。

我想为网络上每个点之间的所有边赋予相同的 ID，这样我就可以计算点之间网络的总长度。我有一个手动解决方案，但这必须使用 > 20,000 点的大型数据集来完成。

install.packages("remotes"); library("remotes")
install_github("luukvdmeer/sfnetworks")
library(sf)
library(sfnetworks)
library(dplyr)

ln <- structure(list(River_ID = c(159, 160, 161, 186, 196), geometry = structure(list(
  structure(c(
    289924.625, 289924.5313, 289922.9688, 289920.0625,
    289915.7499, 289912.7188, 289907.4375, 289905.3438, 289901.1251,
    289889, 289888.5, 289887.5938, 289886.5, 289886.4063, 289885.3124,
    289884.0938, 289884.0001, 289882.8125, 289881.625, 289878.6875,
    289877.9688, 289876.25, 289874.5625, 289874.25, 289872.7188,
    289871.2813, 289871.1875, 289870.0313, 289869, 289868.5939,
    289867.8436, 289865.8438, 289864.0625, 289862.5939, 289862.375,
    289861.5, 289860.7812, 289860.5625, 289859.5313, 289858.375,
    289857.7813, 289855.4063, 289854.25, 289850.8749, 289846.4376,
    289841.9064, 289836.0625, 289828.1562, 289822.8438, 289816.625,
    289812.4376, 289807.9064, 289798.75, 289793.125, 289786.2188,
    289781.375, 289777.3124, 289770.0313, 289765.4375, 289762.2188,
    289759.25, 289755.5938, 289753.0625, 289747.9687, 289743.7499,
    289741.5938, 289739.5, 289736.1874, 289732.75, 289727, 289723.7499,
    289719.625, 289715.5626, 289713.7499, 202817.531300001, 202817.2031,
    202815.1094, 202812.468699999, 202809.3906, 202806.7656,
    202799.7969, 202797.906300001, 202794.093800001, 202783.515699999,
    202783.125, 202782.4844, 202781.906300001, 202781.8125, 202781.3594,
    202781.093800001, 202780.9999, 202780.5469, 202780, 202777.625,
    202777.0469, 202775.718800001, 202774.1875, 202773.906300001,
    202772.1875, 202770.4531, 202770.25, 202768.5156, 202766.6719,
    202766, 202764.0469, 202759.6719, 202755.8749, 202752.781300001,
    202752.1875, 202749.953199999, 202748.297, 202747.906300001,
    202746.0625, 202744.2344, 202743.5625, 202740.4375, 202738.8125,
    202734.5, 202727.9844, 202723.5625, 202719.1875, 202714.9845,
    202713.031300001, 202710.6875, 202710.0469, 202711.406300001,
    202714.5626, 202716.9845, 202718.718900001, 202719.5469,
    202718.734300001, 202716.4531, 202715.125, 202713.7344, 202712.093800001,
    202709.8749, 202708.875, 202709.2655, 202710.7031, 202712.375,
    202712.375, 202712.2344, 202711.0469, 202707.906300001, 202705.406300001,
    202703.0469, 202701.468800001, 202700.7656
  ), .Dim = c(
    74L,
    2L
  ), class = c("XY", "LINESTRING", "sfg")), structure(c(
    289954.375,
    289953.5, 289950.6562, 289949.7499, 289949, 289948.125, 289946.0625,
    289945.9688, 289944.5313, 289943.4063, 289941.3438, 289939.4375,
    289937.4375, 289935.1875, 289932.75, 289930.625, 289928.8125,
    289928.25, 289926.7188, 289925.5313, 289925.7813, 289925.625,
    289925.4063, 289925.1251, 289924.625, 202872.75, 202872.031400001,
    202868.7031, 202867.343699999, 202864.906199999, 202861.515699999,
    202858.297, 202854.406300001, 202851.9375, 202849.468800001,
    202847.703, 202846.75, 202845.4531, 202843.6719, 202843.0625,
    202841.593900001, 202839.7344, 202839.2344, 202838, 202835.9375,
    202832.875, 202825.7344, 202822.9531, 202819.4531, 202817.531300001
  ), .Dim = c(25L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
    290042.6563, 290042.3437, 290041.5313, 290038.4376,
    290037.625, 290036.5313, 290035.5313, 290034.8438, 290034.5313,
    290033.7188, 290032.9375, 290032.125, 290030.3437, 290030.0313,
    290028.625, 290027.5626, 290027.3438, 290026.7188, 290024.5313,
    290023.625, 290020.625, 290018.0001, 290014.9375, 290012.0938,
    290008.5625, 290004.375, 290000.0001, 289999.875, 289997.625,
    289993.7188, 289990.5, 289987.1562, 289985.4063, 289980.375,
    289973.3124, 289966.375, 289961.8438, 289959, 289954.375,
    202884.0625, 202884.25, 202884.843800001, 202888.4531, 202889.75,
    202891.0469, 202892.0469, 202892.656300001, 202892.843800001,
    202893.2501, 202893.5469, 202893.656300001, 202893.4531,
    202893.4531, 202893.343699999, 202893.093800001, 202893.0469,
    202892.843800001, 202891.953199999, 202891.5469, 202889.843800001,
    202888.218800001, 202885.1094, 202880.9219, 202877.5625,
    202873.968800001, 202872.5469, 202872.5156, 202872.625, 202874.5469,
    202876.734300001, 202878.1719, 202877.953199999, 202876.3125,
    202873.468800001, 202872.031400001, 202872.906199999, 202873.0781,
    202872.75
  ), .Dim = c(39L, 2L), class = c(
    "XY", "LINESTRING",
    "sfg"
  )), structure(c(
    290054.125, 290053.4375, 290052.5313,
    290051.625, 290050.0313, 290048.125, 290044.125, 290040.4376,
    290039.4375, 290036.9688, 290031.4375, 290027.5312, 290024.8125,
    290021.7499, 290020.9688, 290018.3437, 290015, 290010.25,
    290006.0313, 290002.4376, 290000.0001, 289999.2187, 289996.6875,
    289995.3438, 289994.125, 289991.1875, 289989.2187, 289987.9688,
    289986.125, 289980.5313, 289975.0314, 289970.9063, 289968.5625,
    289961.0312, 289948.0001, 289939.625, 289933.1563, 289928.3125,
    289926.5313, 289924.625, 202835.953199999, 202835.656300001,
    202835.4531, 202835.343699999, 202835.5469, 202835.7656,
    202836.25, 202836.4531, 202836.5469, 202836.5469, 202835.953199999,
    202836.031400001, 202836.625, 202837.7969, 202838.4844, 202839.343699999,
    202836.25, 202832.7656, 202832.3125, 202833.4844, 202834.4844,
    202834.8125, 202834.2344, 202832.625, 202830.625, 202828.593800001,
    202828.968800001, 202831.0625, 202833.2655, 202835.5781,
    202838, 202838.906199999, 202839.125, 202836.4531, 202830.781300001,
    202827.093800001, 202823.625, 202818.5, 202817.5625, 202817.531300001
  ), .Dim = c(40L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
    290042.625, 290042.0313, 290041.2187, 290040.3125,
    290038.4063, 290037.7188, 290035.8125, 290033.7188, 290030.9063,
    290028.2187, 290021.5313, 290021.2187, 290014.2188, 290013.4063,
    290012.3125, 290010.0625, 290007.9375, 290005.9688, 290004.125,
    290000.0001, 289999.4063, 289998.3125, 289997.5312, 289996.8438,
    289993.625, 289993.0314, 289989.7188, 289989.3438, 289987.625,
    289987.2187, 289984.0313, 289978.125, 289977.9375, 289974.3437,
    289972.7188, 289970.9375, 289967.9375, 289965.2187, 289965.1563,
    289962.3437, 289960.5313, 289959.1251, 289959.0314, 289959.3438,
    289959.4375, 289959.4375, 289959.3438, 289959.2187, 289958.9375,
    289958.5313, 289956.125, 289954.375, 202953.781300001, 202952.4844,
    202951.281300001, 202950.0781, 202948.1875, 202947.5781,
    202945.8749, 202944.281300001, 202941.781300001, 202940.1875,
    202936.375, 202936.1875, 202931.968800001, 202931.4844, 202930.875,
    202929.093800001, 202927.1094, 202925.031300001, 202922.734300001,
    202917.2031, 202916.4375, 202915.2031, 202914.5469, 202914.4531,
    202911.4531, 202910.843800001, 202908.0469, 202907.75, 202906.75,
    202906.5469, 202904.843800001, 202901.843800001, 202901.75,
    202900.0469, 202899.156400001, 202898.0469, 202894.656300001,
    202892.0469, 202891.9844, 202889.343699999, 202887.656300001,
    202885.75, 202884.5469, 202883.343699999, 202882.5469, 202881.343699999,
    202880.0469, 202879.343699999, 202877.656300001, 202876.25,
    202874.25, 202872.75
  ), .Dim = c(52L, 2L), class = c(
    "XY",
    "LINESTRING", "sfg"
  ))
), n_empty = 0L, crs = structure(list(
  input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n    BASEGEOGCRS[\"OSGB 1936\",\n        DATUM[\"OSGB 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n        BBOX[49.75,-9.2,61.14,2.88]],\n    ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
  "sfc_LINESTRING",
  "sfc"
), precision = 0, bbox = structure(c(
  xmin = 289713.7499,
  ymin = 202700.7656, xmax = 290054.125, ymax = 202953.781300001
), class = "bbox"))), row.names = c(NA, -5L), class = c(
  "sf",
  "data.frame"
), sf_column = "geometry", agr = structure(c(River_ID = NA_integer_), .Label = c(
  "constant",
  "aggregate", "identity"
), class = "factor"))

pt <- structure(list(lat = c(
  202805.8942, 202836.136, 202872.9487,
  202905.3284
), lng = c(
  289912.0584, 290014.8446, 290001.2364,
  289984.9382
), id = 1:4, geometry = structure(list(structure(c(
  289912.058400425,
  202805.894199679
), class = c("XY", "POINT", "sfg")), structure(c(
  290014.844597566,
  202836.136003318
), class = c("XY", "POINT", "sfg")), structure(c(
  290001.236395958,
  202872.948712436
), class = c("XY", "POINT", "sfg")), structure(c(
  289984.938209474,
  202905.32838227
), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(
  input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n    BASEGEOGCRS[\"OSGB 1936\",\n        DATUM[\"OSGB 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n        BBOX[49.75,-9.2,61.14,2.88]],\n    ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
  "sfc_POINT",
  "sfc"
), precision = 0, bbox = structure(c(
  xmin = 289912.058400425,
  ymin = 202805.894199679, xmax = 290014.844597566, ymax = 202905.32838227
), class = "bbox"))), row.names = c(NA, 4L), class = c(
  "sf",
  "data.frame"
), sf_column = "geometry", agr = structure(c(
  lat = NA_integer_,
  lng = NA_integer_, 
  id = NA_integer_
), .Label = c(
  "constant",
  "aggregate", "identity"
), class = "factor"))

使它成为一个 sfnetwork

ln_sfnetwork <- as_sfnetwork(ln)

添加点几何

ln_sfnetwork <- st_network_blend(ln_sfnetwork, st_geometry(pt))

通过本质上拆分 sfnetwork 使用点几何创建新边

ln_sf <- ln_sfnetwork %>% 
  activate("edges") %>% 
  mutate(new_river_id = as.character(1:n())) %>% 
  st_as_sf()

用于对节点之间的所有边进行分组的手动位pt。这将出现在点之间的所有边缘分组。

ln_sf[["new_river_id"]][2] <- "1"
ln_sf[["new_river_id"]][c(1, 3, 7, 5, 9)] <- "2"
ln_sf[["new_river_id"]][6] <- "3"
ln_sf[["new_river_id"]][4] <- "4"
ln_sf[["new_river_id"]][8] <- "5"

期望的输出

所需的输出，其中 new_river_id 对于点之间的所有边都相同

然后group_by得到组合长度

ln_sf_merged <- ln_sf %>% 
  group_by(new_river_id) %>% 
  summarise(fraglen = sum(seglen)) %>%
  as.data.frame() %>% 
  select(-geometry)

最后

out_sf <- ln_sf %>%
  left_join(ln_sf_merged, by = "new_river_id")
4

1 回答 1

3

有趣的问题，很高兴看到sfnetworks在道路网络以外的其他类型的网络上使用！我考虑了一下。请注意，答案很长，有时可能很复杂。

很清楚您希望示例的输出是什么，但该示例并未涵盖在混合给定点后网络结构的许多不同可能性。因此，很难概括应该用于合并的规则。但我对此做了一些假设：

据我了解，您希望根据您在河网下游行驶时经过的第一点对边进行分组。因此，如果从边 1 和边 2 向下游（朝向树根）的路径在通过任何其他点之前都通过点 1，则它们应该在同一组中。下游行进到达根而不通过任何点的所有边也应该一起在一个组中。

所以最终目标是为每条边分配一个组索引。在tidygraph（sfnetworks构建的库）中，有几个分组函数可以应用于网络的节点或边缘。此类分组函数旨在用于 tidyverse 动词中，例如mutate()and filter()，其中正在处理的网络是已知的，不需要作为分组函数的参数。因此，您可以运行network %>% mutate(group = group_components())，而无需将网络作为参数显式转发到group_components()。

一切都很好，但据我所知tidygraph，没有实现解决您问题的分组功能。这意味着，我们应该创建自己的！我试了一下。它可能不完全是您所需要的，但至少它应该为您提供一个良好的起点，以便进一步发展。

这个想法是，我们首先计算与您添加到网络中的点相对应的节点到网络中所有节点的距离。然后，我们为每个点选择一组节点，使得从该点到这些节点的旅行时间低于从任何其他点的旅行时间。连接这样一个集合中的节点的边形成属于该点的边组。

我创建了另一个示例网络，它涵盖了更多可能的情况。下面代码的第一部分主要是构建该网络。然后我们开始编写自定义组函数，最后将它应用到我的示例网络以及您的示例网络。

library(sfnetworks)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.0, PROJ 7.2.0
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter

# Create an example network.
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)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

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_2 = sfnetwork(nodes, edges)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid

G = st_network_join(G_1, G_2) %>%
  convert(to_spatial_explicit, .clean = TRUE)

# Create a set of points.
p1 = st_sfc(st_point(c(1, 0.5)))
p2 = st_sfc(st_point(c(2.5, 1.5)))
p3 = st_sfc(st_point(c(1, 2.5)))
p4 = st_sfc(st_point(c(1, 3.5)))
p5 = st_sfc(st_point(c(1.5, 3.5)))
p6 = st_sfc(st_point(c(-1.5, 1.5)))
p7 = st_sfc(st_point(c(-1.5, 2.5)))
p8 = st_sfc(st_point(c(-1.5, 3.5)))
p9 = st_sfc(st_point(c(-1, 0.5)))

# Important!
# Add a column to the points with TRUE values.
# Such that we know which nodes correspond to the points after we inlcude them in the network.
points = st_as_sf(c(p1, p2, p3, p4, p5, p6, p7, p8, p9))
points$is_point = TRUE

plot(G)
plot(st_geometry(points), col = "red", pch = 20, add = TRUE)


# Blend the points into the network.
# Update the is_point column so that all other nodes get a value of FALSE.
G = st_network_blend(G, points) %>%
  morph(to_subgraph, is.na(is_point)) %>%
  mutate(is_point = FALSE) %>%
  unmorph()
#> Warning: st_network_blend assumes attributes are constant over geometries
#> Subsetting by nodes

G
#> # A sfnetwork with 32 nodes and 31 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data:     32 x 2 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -4 ymin: 0 xmax: 4 ymax: 4
#>   geometry is_point
#>    <POINT> <lgl>   
#> 1  (1 0.5) TRUE    
#> 2    (2 1) FALSE   
#> 3    (0 0) FALSE   
#> 4  (1 2.5) TRUE    
#> 5    (1 3) FALSE   
#> 6    (1 2) FALSE   
#> # … with 26 more rows
#> #
#> # Edge Data:     31 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -4 ymin: 0 xmax: 4 ymax: 4
#>    from    to     geometry
#>   <int> <int> <LINESTRING>
#> 1     1     2 (1 0.5, 2 1)
#> 2     3     1 (0 0, 1 0.5)
#> 3     4     5 (1 2.5, 1 3)
#> # … with 28 more rows

plot(G)


# Define our own custom grouping function.
# Remember that in the tidygraph style these functions are meant to be
# used inside tidyverse verbs like mutate() and filter(), in which the
# graph that is being worked on is known and not needed as an input to
# the function.
# The graph being worked on can be accessed with .G().
# Its nodes and edges respectively with .N() and .E().
# The mode argument specifies how to route upstreams on the river network.
# If your network is directed upstreams (i.e. starting at the root of the tree):
# --> Use only outbound edges.
# If your network is directed downstreams (i.e. ending at the root of the tree):
# --> Use only inbound edges.
group_custom = function(mode = "out") {
  # First we get the node indices of the nodes we want to route from.
  # These are:
  # --> All points that were blended into the network.
  # --> The root node at the start of the network tree.
  # Including the root will group all edges that dont have a blended point downstreams.
  origins = which(.N()$is_point | with_graph(.G(), node_is_root()))
  # Calculate the cost matrix from the origins, to all nodes in the network.
  costs = st_network_cost(.G(), from = origins, Inf_as_NaN = TRUE, mode = mode)
  # For each node in the network:
  # --> Define which of the origins is the first to be reached when travelling downstreams.
  # Remember that in the cost matrix:
  # --> The origins (the blended points + the root node) are the rows.
  # --> The destinations (all nodes in the network) are the columns.
  # Hence, we loop over the columns and keep only the minimum cost value per column.
  # We should first remove the zeros, which are the cost values from and to the same node.
  keep_minimum_cost = function(i) {
    i[i == 0] = NaN
    if (any(!is.na(i))) i[i != min(i, na.rm = TRUE)] = NaN
    i
  }
  costs = apply(costs, 2, keep_minimum_cost)
  # For each origin we know now which nodes are in its group.
  # However, we want to know which edges are in the group.
  # The cost matrix does not provide that information.
  # Finding the paths from the origins to the nodes in its group will give us this.
  # Hence, for each origin:
  # --> We compute the paths to all nodes in its group.
  # --> We extract the edge indices that are part of these paths.
  get_edge_indices = function(i) {
    orig = origins[i]
    dest = which(!is.na(costs[i, ]))
    if (length(dest) > 0) {
      paths = st_network_paths(.G(), from = orig, to = dest, mode = mode)
      edge_idxs = do.call(c, paths$edge_paths)
      unique(edge_idxs)
    }
  }
  groups = lapply(seq_len(nrow(costs)), get_edge_indices)
  # In tidygraph the largest group always gets group index number 1.
  # We can achieve the same by ordering the groups by number of edges.
  groups = groups[order(lengths(groups), decreasing = TRUE)]
  # Now we can assign a group index to each edge.
  edge_idxs = do.call(c, groups)
  group_idxs = rep(seq_along(groups), lengths(groups))
  # The last thing left to do is to return the group indices in the correct order.
  # That is: the order of the edges in the edge table of the network.
  group_idxs[order(edge_idxs)]
}

# Group the edges with our custom grouping function.
G = G %>%
  activate("edges") %>%
  mutate(group = group_custom())

G
#> # A sfnetwork with 32 nodes and 31 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Edge Data:     31 x 4 (active)
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: -4 ymin: 0 xmax: 4 ymax: 4
#>    from    to     geometry group
#>   <int> <int> <LINESTRING> <int>
#> 1     1     2 (1 0.5, 2 1)     2
#> 2     3     1 (0 0, 1 0.5)     6
#> 3     4     5 (1 2.5, 1 3)     5
#> 4     6     4 (1 2, 1 2.5)     2
#> 5     6     7   (1 2, 2 3)     2
#> 6     8     9 (1 3.5, 1 4)     7
#> # … with 25 more rows
#> #
#> # Node Data:     32 x 2
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: -4 ymin: 0 xmax: 4 ymax: 4
#>   geometry is_point
#>    <POINT> <lgl>   
#> 1  (1 0.5) TRUE    
#> 2    (2 1) FALSE   
#> 3    (0 0) FALSE   
#> # … with 29 more rows

# Plot the results.
nodes = st_as_sf(G, "nodes")
edges = st_as_sf(G, "edges")
edges$group = as.factor(edges$group) # Such that sf uses categorical colors.

plot(st_geometry(edges))
plot(edges["group"], lwd = 4, key.pos = NULL, add = TRUE)
plot(nodes[nodes$is_point, ], pch = 8, add = TRUE)
plot(nodes[!nodes$is_point, ], pch = 20, add = TRUE)


# Reproduce your example.
# Remember to add the is_point column to the created points.
ln <- structure(list(River_ID = c(159, 160, 161, 186, 196), geometry = structure(list(
  structure(c(
    289924.625, 289924.5313, 289922.9688, 289920.0625,
    289915.7499, 289912.7188, 289907.4375, 289905.3438, 289901.1251,
    289889, 289888.5, 289887.5938, 289886.5, 289886.4063, 289885.3124,
    289884.0938, 289884.0001, 289882.8125, 289881.625, 289878.6875,
    289877.9688, 289876.25, 289874.5625, 289874.25, 289872.7188,
    289871.2813, 289871.1875, 289870.0313, 289869, 289868.5939,
    289867.8436, 289865.8438, 289864.0625, 289862.5939, 289862.375,
    289861.5, 289860.7812, 289860.5625, 289859.5313, 289858.375,
    289857.7813, 289855.4063, 289854.25, 289850.8749, 289846.4376,
    289841.9064, 289836.0625, 289828.1562, 289822.8438, 289816.625,
    289812.4376, 289807.9064, 289798.75, 289793.125, 289786.2188,
    289781.375, 289777.3124, 289770.0313, 289765.4375, 289762.2188,
    289759.25, 289755.5938, 289753.0625, 289747.9687, 289743.7499,
    289741.5938, 289739.5, 289736.1874, 289732.75, 289727, 289723.7499,
    289719.625, 289715.5626, 289713.7499, 202817.531300001, 202817.2031,
    202815.1094, 202812.468699999, 202809.3906, 202806.7656,
    202799.7969, 202797.906300001, 202794.093800001, 202783.515699999,
    202783.125, 202782.4844, 202781.906300001, 202781.8125, 202781.3594,
    202781.093800001, 202780.9999, 202780.5469, 202780, 202777.625,
    202777.0469, 202775.718800001, 202774.1875, 202773.906300001,
    202772.1875, 202770.4531, 202770.25, 202768.5156, 202766.6719,
    202766, 202764.0469, 202759.6719, 202755.8749, 202752.781300001,
    202752.1875, 202749.953199999, 202748.297, 202747.906300001,
    202746.0625, 202744.2344, 202743.5625, 202740.4375, 202738.8125,
    202734.5, 202727.9844, 202723.5625, 202719.1875, 202714.9845,
    202713.031300001, 202710.6875, 202710.0469, 202711.406300001,
    202714.5626, 202716.9845, 202718.718900001, 202719.5469,
    202718.734300001, 202716.4531, 202715.125, 202713.7344, 202712.093800001,
    202709.8749, 202708.875, 202709.2655, 202710.7031, 202712.375,
    202712.375, 202712.2344, 202711.0469, 202707.906300001, 202705.406300001,
    202703.0469, 202701.468800001, 202700.7656
  ), .Dim = c(
    74L,
    2L
  ), class = c("XY", "LINESTRING", "sfg")), structure(c(
    289954.375,
    289953.5, 289950.6562, 289949.7499, 289949, 289948.125, 289946.0625,
    289945.9688, 289944.5313, 289943.4063, 289941.3438, 289939.4375,
    289937.4375, 289935.1875, 289932.75, 289930.625, 289928.8125,
    289928.25, 289926.7188, 289925.5313, 289925.7813, 289925.625,
    289925.4063, 289925.1251, 289924.625, 202872.75, 202872.031400001,
    202868.7031, 202867.343699999, 202864.906199999, 202861.515699999,
    202858.297, 202854.406300001, 202851.9375, 202849.468800001,
    202847.703, 202846.75, 202845.4531, 202843.6719, 202843.0625,
    202841.593900001, 202839.7344, 202839.2344, 202838, 202835.9375,
    202832.875, 202825.7344, 202822.9531, 202819.4531, 202817.531300001
  ), .Dim = c(25L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
    290042.6563, 290042.3437, 290041.5313, 290038.4376,
    290037.625, 290036.5313, 290035.5313, 290034.8438, 290034.5313,
    290033.7188, 290032.9375, 290032.125, 290030.3437, 290030.0313,
    290028.625, 290027.5626, 290027.3438, 290026.7188, 290024.5313,
    290023.625, 290020.625, 290018.0001, 290014.9375, 290012.0938,
    290008.5625, 290004.375, 290000.0001, 289999.875, 289997.625,
    289993.7188, 289990.5, 289987.1562, 289985.4063, 289980.375,
    289973.3124, 289966.375, 289961.8438, 289959, 289954.375,
    202884.0625, 202884.25, 202884.843800001, 202888.4531, 202889.75,
    202891.0469, 202892.0469, 202892.656300001, 202892.843800001,
    202893.2501, 202893.5469, 202893.656300001, 202893.4531,
    202893.4531, 202893.343699999, 202893.093800001, 202893.0469,
    202892.843800001, 202891.953199999, 202891.5469, 202889.843800001,
    202888.218800001, 202885.1094, 202880.9219, 202877.5625,
    202873.968800001, 202872.5469, 202872.5156, 202872.625, 202874.5469,
    202876.734300001, 202878.1719, 202877.953199999, 202876.3125,
    202873.468800001, 202872.031400001, 202872.906199999, 202873.0781,
    202872.75
  ), .Dim = c(39L, 2L), class = c(
    "XY", "LINESTRING",
    "sfg"
  )), structure(c(
    290054.125, 290053.4375, 290052.5313,
    290051.625, 290050.0313, 290048.125, 290044.125, 290040.4376,
    290039.4375, 290036.9688, 290031.4375, 290027.5312, 290024.8125,
    290021.7499, 290020.9688, 290018.3437, 290015, 290010.25,
    290006.0313, 290002.4376, 290000.0001, 289999.2187, 289996.6875,
    289995.3438, 289994.125, 289991.1875, 289989.2187, 289987.9688,
    289986.125, 289980.5313, 289975.0314, 289970.9063, 289968.5625,
    289961.0312, 289948.0001, 289939.625, 289933.1563, 289928.3125,
    289926.5313, 289924.625, 202835.953199999, 202835.656300001,
    202835.4531, 202835.343699999, 202835.5469, 202835.7656,
    202836.25, 202836.4531, 202836.5469, 202836.5469, 202835.953199999,
    202836.031400001, 202836.625, 202837.7969, 202838.4844, 202839.343699999,
    202836.25, 202832.7656, 202832.3125, 202833.4844, 202834.4844,
    202834.8125, 202834.2344, 202832.625, 202830.625, 202828.593800001,
    202828.968800001, 202831.0625, 202833.2655, 202835.5781,
    202838, 202838.906199999, 202839.125, 202836.4531, 202830.781300001,
    202827.093800001, 202823.625, 202818.5, 202817.5625, 202817.531300001
  ), .Dim = c(40L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
    290042.625, 290042.0313, 290041.2187, 290040.3125,
    290038.4063, 290037.7188, 290035.8125, 290033.7188, 290030.9063,
    290028.2187, 290021.5313, 290021.2187, 290014.2188, 290013.4063,
    290012.3125, 290010.0625, 290007.9375, 290005.9688, 290004.125,
    290000.0001, 289999.4063, 289998.3125, 289997.5312, 289996.8438,
    289993.625, 289993.0314, 289989.7188, 289989.3438, 289987.625,
    289987.2187, 289984.0313, 289978.125, 289977.9375, 289974.3437,
    289972.7188, 289970.9375, 289967.9375, 289965.2187, 289965.1563,
    289962.3437, 289960.5313, 289959.1251, 289959.0314, 289959.3438,
    289959.4375, 289959.4375, 289959.3438, 289959.2187, 289958.9375,
    289958.5313, 289956.125, 289954.375, 202953.781300001, 202952.4844,
    202951.281300001, 202950.0781, 202948.1875, 202947.5781,
    202945.8749, 202944.281300001, 202941.781300001, 202940.1875,
    202936.375, 202936.1875, 202931.968800001, 202931.4844, 202930.875,
    202929.093800001, 202927.1094, 202925.031300001, 202922.734300001,
    202917.2031, 202916.4375, 202915.2031, 202914.5469, 202914.4531,
    202911.4531, 202910.843800001, 202908.0469, 202907.75, 202906.75,
    202906.5469, 202904.843800001, 202901.843800001, 202901.75,
    202900.0469, 202899.156400001, 202898.0469, 202894.656300001,
    202892.0469, 202891.9844, 202889.343699999, 202887.656300001,
    202885.75, 202884.5469, 202883.343699999, 202882.5469, 202881.343699999,
    202880.0469, 202879.343699999, 202877.656300001, 202876.25,
    202874.25, 202872.75
  ), .Dim = c(52L, 2L), class = c(
    "XY",
    "LINESTRING", "sfg"
  ))
), n_empty = 0L, crs = structure(list(
  input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n    BASEGEOGCRS[\"OSGB 1936\",\n        DATUM[\"OSGB 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n        BBOX[49.75,-9.2,61.14,2.88]],\n    ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
  "sfc_LINESTRING",
  "sfc"
), precision = 0, bbox = structure(c(
  xmin = 289713.7499,
  ymin = 202700.7656, xmax = 290054.125, ymax = 202953.781300001
), class = "bbox"))), row.names = c(NA, -5L), class = c(
  "sf",
  "data.frame"
), sf_column = "geometry", agr = structure(c(River_ID = NA_integer_), .Label = c(
  "constant",
  "aggregate", "identity"
), class = "factor"))

pt <- structure(list(lat = c(
  202805.8942, 202836.136, 202872.9487,
  202905.3284
), lng = c(
  289912.0584, 290014.8446, 290001.2364,
  289984.9382
), id = 1:4, geometry = structure(list(structure(c(
  289912.058400425,
  202805.894199679
), class = c("XY", "POINT", "sfg")), structure(c(
  290014.844597566,
  202836.136003318
), class = c("XY", "POINT", "sfg")), structure(c(
  290001.236395958,
  202872.948712436
), class = c("XY", "POINT", "sfg")), structure(c(
  289984.938209474,
  202905.32838227
), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(
  input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n    BASEGEOGCRS[\"OSGB 1936\",\n        DATUM[\"OSGB 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n        BBOX[49.75,-9.2,61.14,2.88]],\n    ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
  "sfc_POINT",
  "sfc"
), precision = 0, bbox = structure(c(
  xmin = 289912.058400425,
  ymin = 202805.894199679, xmax = 290014.844597566, ymax = 202905.32838227
), class = "bbox"))), row.names = c(NA, 4L), class = c(
  "sf",
  "data.frame"
), sf_column = "geometry", agr = structure(c(
  lat = NA_integer_,
  lng = NA_integer_, 
  id = NA_integer_
), .Label = c(
  "constant",
  "aggregate", "identity"
), class = "factor"))

ln_sfnetwork <- as_sfnetwork(ln)

pt$is_point <- TRUE

ln_sfnetwork <- st_network_blend(ln_sfnetwork, pt["is_point"]) %>%
  morph(to_subgraph, is.na(is_point)) %>%
  mutate(is_point = FALSE) %>%
  unmorph()
#> Warning: st_network_blend assumes attributes are constant over geometries
#> Subsetting by nodes

# NOTE! Your network is directed towards the root of the tree!
# Therefore we route over inbound edges instead of outbound edges to get correct results.
ln_sfnetwork <- ln_sfnetwork %>%
  activate("edges") %>%
  mutate(group = group_custom(mode = "in"))

# Plot the results.
nodes = st_as_sf(ln_sfnetwork, "nodes")
edges = st_as_sf(ln_sfnetwork, "edges")
edges$group = as.factor(edges$group) # Such that sf uses categorical colors.

plot(st_geometry(edges))
plot(edges["group"], lwd = 4, key.pos = NULL, add = TRUE)
plot(nodes[nodes$is_point, ], pch = 8, add = TRUE)
plot(nodes[!nodes$is_point, ], pch = 20, add = TRUE)

由reprex 包（v0.3.0）于 2021 年 2 月 3 日创建

于 2021-02-03T20:41:47.647 回答