问题描述
我有一个有根的树，带有空间显式边 (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.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),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.343699999,202893.093800001,202893.0469,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.906199999,202873.0781,202872.75
),.Dim = c(39L,class = c(
"XY","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,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,202835.953199999,202835.656300001,202835.4531,202835.343699999,202835.5469,202835.7656,202836.25,202836.4531,202836.5469,202836.031400001,202836.625,202837.7969,202838.4844,202839.343699999,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.906199999,202839.125,202830.781300001,202827.093800001,202823.625,202818.5,202817.5625,.Dim = c(40L,structure(c(
290042.625,290042.0313,290041.2187,290040.3125,290038.4063,290037.7188,290035.8125,290030.9063,290028.2187,290021.5313,290021.2187,290014.2188,290013.4063,290012.3125,290010.0625,290007.9375,290005.9688,290004.125,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.2187,289958.9375,289958.5313,289956.125,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,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,.Dim = c(52L,"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\",\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,0.0174532925199433],8801]],\n ParaMETER[\"Longitude of natural origin\",-2,8802]],\n ParaMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],8805]],\n ParaMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",8806]],\n ParaMETER[\"False northing\",-100000,8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],1]],\n AXIS[\"(N)\",north,\n ORDER[2],\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
),"POINT",structure(c(
290014.844597566,202836.136003318
),structure(c(
290001.236395958,202872.948712436
),structure(c(
289984.938209474,202905.32838227
),"sfg"))),class = c(
"sfc_POINT",bBox = structure(c(
xmin = 289912.058400425,ymin = 202805.894199679,xmax = 290014.844597566,ymax = 202905.32838227
),4L),agr = structure(c(
lat = NA_integer_,lng = NA_integer_,id = NA_integer_
),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"
期望的输出
然后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 的库）中，有几个分组函数可以应用于网络的节点或边。此类分组函数旨在用于诸如 mutate() 和 filter() 之类的 tidyverse 动词，其中正在处理的网络是已知的，不需要作为分组函数的参数。因此，您可以运行 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,3,5,8,9,9)
to = c(5,6,4,7,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,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,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,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.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),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.343699999,202893.093800001,202893.0469,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.906199999,202873.0781,202872.75
),.Dim = c(39L,class = c(
"XY","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,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,202835.953199999,202835.656300001,202835.4531,202835.343699999,202835.5469,202835.7656,202836.25,202836.4531,202836.5469,202836.031400001,202836.625,202837.7969,202838.4844,202839.343699999,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.906199999,202839.125,202830.781300001,202827.093800001,202823.625,202818.5,202817.5625,.Dim = c(40L,structure(c(
290042.625,290042.0313,290041.2187,290040.3125,290038.4063,290037.7188,290035.8125,290030.9063,290028.2187,290021.5313,290021.2187,290014.2188,290013.4063,290012.3125,290010.0625,290007.9375,290005.9688,290004.125,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.2187,289958.9375,289958.5313,289956.125,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,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,.Dim = c(52L,"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\",\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,0.0174532925199433],8801]],\n PARAMETER[\"Longitude of natural origin\",-2,8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],8805]],\n PARAMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",8806]],\n PARAMETER[\"False northing\",-100000,8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],1]],\n AXIS[\"(N)\",north,\n ORDER[2],\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
),"POINT",structure(c(
290014.844597566,202836.136003318
),structure(c(
290001.236395958,202872.948712436
),structure(c(
289984.938209474,202905.32838227
),"sfg"))),class = c(
"sfc_POINT",bbox = structure(c(
xmin = 289912.058400425,ymin = 202805.894199679,xmax = 290014.844597566,ymax = 202905.32838227
),4L),agr = structure(c(
lat = NA_integer_,lng = NA_integer_,id = NA_integer_
),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,add = TRUE)
由 reprex package (v0.3.0) 于 2021 年 2 月 3 日创建