使用 r 中的 sfnetwork 点之间所有边的组 ID
Group ID for all edges between points using sfnetwork in r
我有一个带有空间显式边 (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))
使用点 geom 通过本质上分割 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"
期望的输出
然后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")
有趣的问题，很高兴看到 sfnetworks 被用于道路网络以外的其他类型的网络！我考虑了一下。当心，答案很长，有时可能很复杂。
很清楚您希望示例的输出是什么，但是该示例并未涵盖混合给定点后网络结构的许多不同可能性。因此很难概括应该用于合并的规则。但我对此做了一些假设：
据我了解，您想根据 第一个点 是您在下游河网。因此，如果从边 1 和边 2 向下游（朝向树的根）的路径在通过任何其他点之前都经过点 1，则它们应该在同一组中。下游行进到达根而不经过任何点的所有边也应该在一个组中。
所以最终目标是为每条边分配一个组索引。在 tidygraph（构建 sfnetworks 的库）中，有几个分组函数可以应用于网络的节点或边缘。此类分组函数旨在用于 tidyverse 动词，如 mutate() 和 filter()，其中正在处理的网络是已知的，不需要作为分组函数的参数。因此，您可以 运行 network %>% mutate(group = group_components())，而无需显式将网络作为参数转发给 group_components().
一切都很好，但据我所知 tidygraph 没有实施分组功能来解决您的问题。这意味着，我们应该创建自己的！我试了一下。它可能不完全是您所需要的，但至少它应该为您提供一个良好的起点以进一步构建。
思路是我们先计算到你添加到网络的点对应的节点的距离，到所有网络中的节点。然后，我们为每个点 select 一组节点，使得从该点到这些节点的行进时间比从任何其他点的行进时间 低 。连接这样一个集合中的节点的边形成属于该点的边组。
我创建了另一个示例网络，它涵盖了更多可能的情况。下面代码的第一部分主要是构建该网络。然后我们开始编写自定义组函数，最后将其应用到我的示例网络以及您的示例网络。
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 package (v0.3.0)
于 2021-02-03 创建
我有一个带有空间显式边 (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))
使用点 geom 通过本质上分割 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"
期望的输出
然后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")
有趣的问题，很高兴看到 sfnetworks 被用于道路网络以外的其他类型的网络！我考虑了一下。当心，答案很长，有时可能很复杂。
很清楚您希望示例的输出是什么，但是该示例并未涵盖混合给定点后网络结构的许多不同可能性。因此很难概括应该用于合并的规则。但我对此做了一些假设：
据我了解，您想根据 第一个点 是您在下游河网。因此，如果从边 1 和边 2 向下游（朝向树的根）的路径在通过任何其他点之前都经过点 1，则它们应该在同一组中。下游行进到达根而不经过任何点的所有边也应该在一个组中。
所以最终目标是为每条边分配一个组索引。在 tidygraph（构建 sfnetworks 的库）中，有几个分组函数可以应用于网络的节点或边缘。此类分组函数旨在用于 tidyverse 动词，如 mutate() 和 filter()，其中正在处理的网络是已知的，不需要作为分组函数的参数。因此，您可以 运行 network %>% mutate(group = group_components())，而无需显式将网络作为参数转发给 group_components().
一切都很好，但据我所知 tidygraph 没有实施分组功能来解决您的问题。这意味着，我们应该创建自己的！我试了一下。它可能不完全是您所需要的，但至少它应该为您提供一个良好的起点以进一步构建。
思路是我们先计算到你添加到网络的点对应的节点的距离，到所有网络中的节点。然后，我们为每个点 select 一组节点，使得从该点到这些节点的行进时间比从任何其他点的行进时间 低 。连接这样一个集合中的节点的边形成属于该点的边组。
我创建了另一个示例网络，它涵盖了更多可能的情况。下面代码的第一部分主要是构建该网络。然后我们开始编写自定义组函数，最后将其应用到我的示例网络以及您的示例网络。
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 package (v0.3.0)
于 2021-02-03 创建