如何将 sf 空间点列表转换为可路由图
How to convert a list of sf spatial points into a routable graph
我有一个 sf dataframe 对象，其中包含一系列代表公交路线形状的点。我想把这个对象变成一个可路由的图，这样我就可以估计从 c 点遍历到 t.
所需的时间
这是我使用 dodgr package 尝试过的方法，但我不确定我在这里做错了什么：
library(dodgr)
graph <- weight_streetnet(mydata, wt_profile = "motorcar", type_col="highway" , id_col = "id")
Error in check_highway_osmid(x, wt_profile) :
Please specify type_col to be used for weighting streetnet
可重现的数据
数据如下图所示
mydata <- structure(list(shape_id = c(52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L), length = structure(c(0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197), units = structure(list(
numerator = "km", denominator = character(0)), class = "symbolic_units"), class = "units"),
geometry = structure(list(structure(c(-46.5623281998182,
-23.5213458001468), class = c("XY", "POINT", "sfg")), structure(c(-46.562221,
-23.52129), class = c("XY", "POINT", "sfg")), structure(c(-46.562121,
-23.521235), class = c("XY", "POINT", "sfg")), structure(c(-46.5620233332577,
-23.5211840000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561925666591,
-23.5211330000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561828,
-23.521082), class = c("XY", "POINT", "sfg")), structure(c(-46.5618098335317,
-23.5212126666783), class = c("XY", "POINT", "sfg")), structure(c(-46.5617916670273,
-23.5213433333544), class = c("XY", "POINT", "sfg")), structure(c(-46.5617735004869,
-23.5214740000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5617553339104,
-23.5216046667004), class = c("XY", "POINT", "sfg")), structure(c(-46.5617371672978,
-23.5217353333702), class = c("XY", "POINT", "sfg")), structure(c(-46.5617190006492,
-23.5218660000379), class = c("XY", "POINT", "sfg")), structure(c(-46.5617008339645,
-23.5219966667036), class = c("XY", "POINT", "sfg")), structure(c(-46.5616826672438,
-23.5221273333671), class = c("XY", "POINT", "sfg")), structure(c(-46.5616645004869,
-23.5222580000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5616463336941,
-23.5223886666877), class = c("XY", "POINT", "sfg")), structure(c(-46.5616281668651,
-23.5225193333449), class = c("XY", "POINT", "sfg")), structure(c(-46.56161,
-23.52265), class = c("XY", "POINT", "sfg")), structure(c(-46.5617355000207,
-23.5226427501509), class = c("XY", "POINT", "sfg")), structure(c(-46.5618610000276,
-23.5226355002012), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = -46.5623281998182,
ymin = -23.52265, xmax = -46.56161, ymax = -23.521082), class = "bbox"), crs = structure(list(
epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L),
id = c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j",
"k", "l", "m", "n", "o", "p", "q", "r", "s", "t"), speed_kmh = c(11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 11)), sf_column = "geometry", agr = structure(c(shape_id = NA_integer_,
length = NA_integer_, id = NA_integer_, speed_kmh = NA_integer_
), class = "factor", .Label = c("constant", "aggregate", "identity"
)), row.names = c("1.13", "1.14", "1.15", "1.16", "1.17", "1.18",
"1.19", "1.20", "1.21", "1.22", "1.23", "1.24", "1.25", "1.26",
"1.27", "1.28", "1.29", "1.30", "1.31", "1.32"), class = c("sf",
"data.table", "data.frame"))
我认为您可以通过将数据转换为 igraph 对象并使用 igraph 库中的功能来解决它。
您需要建立 Edges 和 Vertex 以及 weight 值。
在 igraph 中，边是 link 表示两个节点（源和目标）之间的连接。在这种情况下，link 是 "street"，点是节点。
library(igraph)
GraphResult <- data.frame(Source = c(NULL),
Target = c(NULL),
weight = c(NULL))
for (i in 1:(dim(mydata)[1] - 1)) {
TempGraphResult <- data.frame(Source = c(0),
Target = c(0),
weight = c(0))
TempGraphResult$Source[1] <- mydata$id[i]
TempGraphResult$Target[1] <- mydata$id[i + 1]
TempGraphResult$weight[1] <- mydata$length[i]
GraphResult <- rbind(GraphResult, TempGraphResult) }
MyIgraph <- graph_from_data_frame(GraphResult)
#In this case works perfectly. But if you have more weight variables and even
#additional variables for the nodes, igraph have functions for constructing the
#igraph object
distances(MyIgraph, "c", "t") #returns 3.254183. Seems correct (0.1914225*17)
SquareMatrix <- distances(MyIgraph)
#*distances() is a function from igraph that performs the routing calculations.
可以实现更复杂的网络和计算路由。例如，您可以设置道路的方向。
也许 dodger 可以解决这个问题，但我不确定。
weight_streetnet 函数实际上仅设计用于处理实际街道网络，通常由 osmdata::osmdata_sf/sp/sc() 函数生成。尽管如此，还是可以对其进行调整以处理此类情况。需要做的主要事情是将点转换成知道它们之间的边的东西，比如 sf::LINESTRING 对象：
x <- sf::st_combine (mydata) %>%
sf::st_cast ("LINESTRING") %>%
sf::st_sf ()
这给出了一个单行对象，然后可以将其转换为 dodgr 格式，并且 id 值匹配回边缘
net <- weight_streetnet (x, type_col = "shape_id", id_col = "id", wt_profile = 1)
net$from_id <- mydata$id [as.integer (net$from_id)]
net$to_id <- mydata$id [as.integer (net$to_id)]
到那时，dodgr 将直接根据地理坐标计算并插入距离。您的距离也可以通过替换 d_weighted 值插入并用于路由：
net$d_weighted <- as.numeric (mydata$length [1])
dodgr_dists (net, from = "c", to = "t") # 236.0481
如果您真的希望您的距离代表用于计算最终结果的绝对距离，那么只需替换 $d 值
net$d <- net$d_weighted
dodgr_dists (net, from = "c", to = "t") # 3.254183
请注意，对于 "simple" 这样的问题，igraph 通常会更快，因为它使用一组权重计算路线。在这种情况下 dodgr 的唯一真正优势是能够使用 "dual weights" - $d_weighted 和 $d 值 - 这样根据 [=24= 计算路线], 最终距离根据 $d.
如果您想将其包含在 'tidy' 工作流程中，您还可以考虑混合使用 sf 和 tidygraph。后者以 tbl_graph class 的形式为 networks/graphs 提供了一个简洁的框架，其中子 class 是 igraph（因此，您可以使用 [=所有 igraph 中的 22=] 个对象作为一个 igraph 对象）。但是，您可以将节点和边分析为 tibbles，并使用 filter()、select()、mutate() 等函数。当然，这些tibbles也可以包含我们从sf知道的几何列表列，为节点和边添加地理信息。
该方法远非完美，欢迎进行改进，但它仍然展示了另一种处理问题的方法。
# Load libraries.
library(tidyverse)
library(sf)
library(tidygraph)
library(igraph)
library(units)
就像其他答案一样，我们需要在节点之间创建边。现在，我假设这些点只是按字母顺序连接。但是，对于 tidygraph 方法，我们似乎需要数字 ID 而不是字符。
# Add a numeric ID column to the nodes.
nodes <- mydata %>%
rename(id_chr = id) %>%
rowid_to_column("id") %>%
select(id, id_chr, everything())
# Define the source node of each edge, and the target node of each edge.
sources <- nodes %>% slice(-n())
targets <- nodes %>% slice(-1)
# Write a function to create lines between data frames of source and target points.
pt2l <- function(x, y) { st_linestring(rbind(st_coordinates(x), st_coordinates(y))) }
# Create the edges.
edges <- tibble(
from = sources %>% pull(id),
to = targets %>% pull(id),
length = sources %>% pull(length),
speed = sources %>% pull(speed_kmh),
geometry = map2(st_geometry(sources), st_geometry(targets), pt2l)
) %>% st_as_sf() %>% st_set_crs(st_crs(nodes))
# Add a time column to the edges.
edges <- edges %>%
mutate(speed = set_units(speed, "km/h")) %>%
mutate(time = length / speed)
# Clean up the nodes data.
nodes <- nodes %>%
select(-length, -speed_kmh)
# Create the tbl_graph object out of the nodes and edges.
# Providing the edges as sf object is problematic for tidygraph, unfortunately.
# Therefore, we have to provide them as a tibble.
graph <- tbl_graph(nodes = nodes, edges = as_tibble(edges), directed = FALSE)
这为我们提供了以下 tbl_graph 对象：
# A tbl_graph: 20 nodes and 19 edges
#
# An undirected simple graph with 1 component
#
# Node Data: 20 x 4 (active)
id id_chr shape_id geometry
<int> <chr> <int> <POINT [°]>
1 1 a 52421 (-46.56233 -23.52135)
2 2 b 52421 (-46.56222 -23.52129)
3 3 c 52421 (-46.56212 -23.52124)
4 4 d 52421 (-46.56202 -23.52118)
5 5 e 52421 (-46.56193 -23.52113)
6 6 f 52421 (-46.56183 -23.52108)
# … with 14 more rows
#
# Edge Data: 19 x 6
from to length speed geometry time
<int> <int> [km] [km/h] <LINESTRING [°]> [h]
1 1 2 0.1914225 11 (-46.56233 -23.52135, -46.56222 -23.5… 0.017402…
2 2 3 0.1914225 11 (-46.56222 -23.52129, -46.56212 -23.5… 0.017402…
3 3 4 0.1914225 11 (-46.56212 -23.52124, -46.56202 -23.5… 0.017402…
# … with 16 more rows
现在我们在图形结构中拥有所有内容，我们可以 select 我们想要路由的节点和我们想要路由到的节点，并找到它们之间的最短路径，行程时间为权重变量，使用 igraph 中的 shortest_path 函数。我们现在只使用一对一路由（'c' 到 't'），但对于一对多、多对一或多对-很多。
# Select the node from which and to which the shortest path should be found.
from_node <- graph %>%
activate(nodes) %>%
filter(id_chr == "c") %>%
pull(id)
to_node <- graph %>%
activate(nodes) %>%
filter(id_chr == "t") %>%
pull(id)
# Find the shortest path between these nodes
path <- shortest_paths(
graph = graph,
from = from_node,
to = to_node,
output = 'both',
weights = graph %>% activate(edges) %>% pull(time)
)
生成的路径是一个列表，其中包含构成路径的节点和边。
$vpath
$vpath[[1]]
+ 18/20 vertices, from e43a089:
[1] 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
$epath
$epath[[1]]
+ 17/19 edges from e43a089:
[1] 3-- 4 4-- 5 5-- 6 6-- 7 7-- 8 8-- 9 9--10 10--11 11--12 12--13
[11] 13--14 14--15 15--16 16--17 17--18 18--19 19--20
我们可以创建一个只包含最短路径的节点和边的原始图的子图。
path_graph <- graph %>%
subgraph.edges(eids = path$epath %>% unlist()) %>%
as_tbl_graph()
# A tbl_graph: 18 nodes and 17 edges
#
# An undirected simple graph with 1 component
#
# Node Data: 18 x 4 (active)
id id_chr shape_id geometry
<int> <chr> <int> <POINT [°]>
1 3 c 52421 (-46.56212 -23.52124)
2 4 d 52421 (-46.56202 -23.52118)
3 5 e 52421 (-46.56193 -23.52113)
4 6 f 52421 (-46.56183 -23.52108)
5 7 g 52421 (-46.56181 -23.52121)
6 8 h 52421 (-46.56179 -23.52134)
# … with 12 more rows
#
# Edge Data: 17 x 6
from to length speed geometry time
<int> <int> [km] [km/h] <LINESTRING [°]> [h]
1 1 2 0.1914225 11 (-46.56212 -23.52124, -46.56202 -23.5… 0.017402…
2 2 3 0.1914225 11 (-46.56202 -23.52118, -46.56193 -23.5… 0.017402…
3 3 4 0.1914225 11 (-46.56193 -23.52113, -46.56183 -23.5… 0.017402…
# … with 14 more rows
在这里，发生了一些我不喜欢的事情。 Tidygraph/igraph 似乎有一个内部节点 ID 结构，你看到在子图中，egdes 数据中的 from 和 to 列与我们的 id 不匹配不再是节点数据中的列，而是简单地引用节点数据的行号。我不确定如何解决这个问题。
无论哪种方式，我们现在都有从 'c' 到 't' 的路径作为子图，并且可以轻松地对其进行分析。例如，通过计算路径的总行程时间（就像问题一样）。
path_graph %>%
activate(edges) %>%
as_tibble() %>%
summarise(total_time = sum(time))
# A tibble: 1 x 1
total_time
[h]
1 0.2958348
但绘制它也很容易，保留了地理信息（只需将节点和边导出为 sf 对象）。
ggplot() +
geom_sf(data = graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), col = 'darkgrey') +
geom_sf(data = graph %>% activate(nodes) %>% as_tibble() %>% st_as_sf(), col = 'darkgrey', size = 0.5) +
geom_sf(data = path_graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), lwd = 1, col = 'firebrick') +
geom_sf(data = path_graph %>% activate(nodes) %>% filter(id %in% c(from_node, to_node)) %>% as_tibble() %>% st_as_sf(), size = 2)
plot
关于这种 tidygraph-sf 方法的 r-spatial 博客 post 可能会出现 ;)
我有一个 sf dataframe 对象，其中包含一系列代表公交路线形状的点。我想把这个对象变成一个可路由的图，这样我就可以估计从 c 点遍历到 t.
这是我使用 dodgr package 尝试过的方法，但我不确定我在这里做错了什么：
library(dodgr)
graph <- weight_streetnet(mydata, wt_profile = "motorcar", type_col="highway" , id_col = "id")
Error in check_highway_osmid(x, wt_profile) : Please specify type_col to be used for weighting streetnet
可重现的数据
数据如下图所示
mydata <- structure(list(shape_id = c(52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L), length = structure(c(0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197), units = structure(list(
numerator = "km", denominator = character(0)), class = "symbolic_units"), class = "units"),
geometry = structure(list(structure(c(-46.5623281998182,
-23.5213458001468), class = c("XY", "POINT", "sfg")), structure(c(-46.562221,
-23.52129), class = c("XY", "POINT", "sfg")), structure(c(-46.562121,
-23.521235), class = c("XY", "POINT", "sfg")), structure(c(-46.5620233332577,
-23.5211840000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561925666591,
-23.5211330000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561828,
-23.521082), class = c("XY", "POINT", "sfg")), structure(c(-46.5618098335317,
-23.5212126666783), class = c("XY", "POINT", "sfg")), structure(c(-46.5617916670273,
-23.5213433333544), class = c("XY", "POINT", "sfg")), structure(c(-46.5617735004869,
-23.5214740000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5617553339104,
-23.5216046667004), class = c("XY", "POINT", "sfg")), structure(c(-46.5617371672978,
-23.5217353333702), class = c("XY", "POINT", "sfg")), structure(c(-46.5617190006492,
-23.5218660000379), class = c("XY", "POINT", "sfg")), structure(c(-46.5617008339645,
-23.5219966667036), class = c("XY", "POINT", "sfg")), structure(c(-46.5616826672438,
-23.5221273333671), class = c("XY", "POINT", "sfg")), structure(c(-46.5616645004869,
-23.5222580000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5616463336941,
-23.5223886666877), class = c("XY", "POINT", "sfg")), structure(c(-46.5616281668651,
-23.5225193333449), class = c("XY", "POINT", "sfg")), structure(c(-46.56161,
-23.52265), class = c("XY", "POINT", "sfg")), structure(c(-46.5617355000207,
-23.5226427501509), class = c("XY", "POINT", "sfg")), structure(c(-46.5618610000276,
-23.5226355002012), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = -46.5623281998182,
ymin = -23.52265, xmax = -46.56161, ymax = -23.521082), class = "bbox"), crs = structure(list(
epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L),
id = c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j",
"k", "l", "m", "n", "o", "p", "q", "r", "s", "t"), speed_kmh = c(11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 11)), sf_column = "geometry", agr = structure(c(shape_id = NA_integer_,
length = NA_integer_, id = NA_integer_, speed_kmh = NA_integer_
), class = "factor", .Label = c("constant", "aggregate", "identity"
)), row.names = c("1.13", "1.14", "1.15", "1.16", "1.17", "1.18",
"1.19", "1.20", "1.21", "1.22", "1.23", "1.24", "1.25", "1.26",
"1.27", "1.28", "1.29", "1.30", "1.31", "1.32"), class = c("sf",
"data.table", "data.frame"))
我认为您可以通过将数据转换为 igraph 对象并使用 igraph 库中的功能来解决它。 您需要建立 Edges 和 Vertex 以及 weight 值。 在 igraph 中，边是 link 表示两个节点（源和目标）之间的连接。在这种情况下，link 是 "street"，点是节点。
library(igraph)
GraphResult <- data.frame(Source = c(NULL),
Target = c(NULL),
weight = c(NULL))
for (i in 1:(dim(mydata)[1] - 1)) {
TempGraphResult <- data.frame(Source = c(0),
Target = c(0),
weight = c(0))
TempGraphResult$Source[1] <- mydata$id[i]
TempGraphResult$Target[1] <- mydata$id[i + 1]
TempGraphResult$weight[1] <- mydata$length[i]
GraphResult <- rbind(GraphResult, TempGraphResult) }
MyIgraph <- graph_from_data_frame(GraphResult)
#In this case works perfectly. But if you have more weight variables and even
#additional variables for the nodes, igraph have functions for constructing the
#igraph object
distances(MyIgraph, "c", "t") #returns 3.254183. Seems correct (0.1914225*17)
SquareMatrix <- distances(MyIgraph)
#*distances() is a function from igraph that performs the routing calculations.
可以实现更复杂的网络和计算路由。例如，您可以设置道路的方向。
也许 dodger 可以解决这个问题，但我不确定。
weight_streetnet 函数实际上仅设计用于处理实际街道网络，通常由 osmdata::osmdata_sf/sp/sc() 函数生成。尽管如此，还是可以对其进行调整以处理此类情况。需要做的主要事情是将点转换成知道它们之间的边的东西，比如 sf::LINESTRING 对象：
x <- sf::st_combine (mydata) %>%
sf::st_cast ("LINESTRING") %>%
sf::st_sf ()
这给出了一个单行对象，然后可以将其转换为 dodgr 格式，并且 id 值匹配回边缘
net <- weight_streetnet (x, type_col = "shape_id", id_col = "id", wt_profile = 1)
net$from_id <- mydata$id [as.integer (net$from_id)]
net$to_id <- mydata$id [as.integer (net$to_id)]
到那时，dodgr 将直接根据地理坐标计算并插入距离。您的距离也可以通过替换 d_weighted 值插入并用于路由：
net$d_weighted <- as.numeric (mydata$length [1])
dodgr_dists (net, from = "c", to = "t") # 236.0481
如果您真的希望您的距离代表用于计算最终结果的绝对距离，那么只需替换 $d 值
net$d <- net$d_weighted
dodgr_dists (net, from = "c", to = "t") # 3.254183
请注意，对于 "simple" 这样的问题，igraph 通常会更快，因为它使用一组权重计算路线。在这种情况下 dodgr 的唯一真正优势是能够使用 "dual weights" - $d_weighted 和 $d 值 - 这样根据 [=24= 计算路线], 最终距离根据 $d.
如果您想将其包含在 'tidy' 工作流程中，您还可以考虑混合使用 sf 和 tidygraph。后者以 tbl_graph class 的形式为 networks/graphs 提供了一个简洁的框架，其中子 class 是 igraph（因此，您可以使用 [=所有 igraph 中的 22=] 个对象作为一个 igraph 对象）。但是，您可以将节点和边分析为 tibbles，并使用 filter()、select()、mutate() 等函数。当然，这些tibbles也可以包含我们从sf知道的几何列表列，为节点和边添加地理信息。
该方法远非完美，欢迎进行改进，但它仍然展示了另一种处理问题的方法。
# Load libraries.
library(tidyverse)
library(sf)
library(tidygraph)
library(igraph)
library(units)
就像其他答案一样，我们需要在节点之间创建边。现在，我假设这些点只是按字母顺序连接。但是，对于 tidygraph 方法，我们似乎需要数字 ID 而不是字符。
# Add a numeric ID column to the nodes.
nodes <- mydata %>%
rename(id_chr = id) %>%
rowid_to_column("id") %>%
select(id, id_chr, everything())
# Define the source node of each edge, and the target node of each edge.
sources <- nodes %>% slice(-n())
targets <- nodes %>% slice(-1)
# Write a function to create lines between data frames of source and target points.
pt2l <- function(x, y) { st_linestring(rbind(st_coordinates(x), st_coordinates(y))) }
# Create the edges.
edges <- tibble(
from = sources %>% pull(id),
to = targets %>% pull(id),
length = sources %>% pull(length),
speed = sources %>% pull(speed_kmh),
geometry = map2(st_geometry(sources), st_geometry(targets), pt2l)
) %>% st_as_sf() %>% st_set_crs(st_crs(nodes))
# Add a time column to the edges.
edges <- edges %>%
mutate(speed = set_units(speed, "km/h")) %>%
mutate(time = length / speed)
# Clean up the nodes data.
nodes <- nodes %>%
select(-length, -speed_kmh)
# Create the tbl_graph object out of the nodes and edges.
# Providing the edges as sf object is problematic for tidygraph, unfortunately.
# Therefore, we have to provide them as a tibble.
graph <- tbl_graph(nodes = nodes, edges = as_tibble(edges), directed = FALSE)
这为我们提供了以下 tbl_graph 对象：
# A tbl_graph: 20 nodes and 19 edges
#
# An undirected simple graph with 1 component
#
# Node Data: 20 x 4 (active)
id id_chr shape_id geometry
<int> <chr> <int> <POINT [°]>
1 1 a 52421 (-46.56233 -23.52135)
2 2 b 52421 (-46.56222 -23.52129)
3 3 c 52421 (-46.56212 -23.52124)
4 4 d 52421 (-46.56202 -23.52118)
5 5 e 52421 (-46.56193 -23.52113)
6 6 f 52421 (-46.56183 -23.52108)
# … with 14 more rows
#
# Edge Data: 19 x 6
from to length speed geometry time
<int> <int> [km] [km/h] <LINESTRING [°]> [h]
1 1 2 0.1914225 11 (-46.56233 -23.52135, -46.56222 -23.5… 0.017402…
2 2 3 0.1914225 11 (-46.56222 -23.52129, -46.56212 -23.5… 0.017402…
3 3 4 0.1914225 11 (-46.56212 -23.52124, -46.56202 -23.5… 0.017402…
# … with 16 more rows
现在我们在图形结构中拥有所有内容，我们可以 select 我们想要路由的节点和我们想要路由到的节点，并找到它们之间的最短路径，行程时间为权重变量，使用 igraph 中的 shortest_path 函数。我们现在只使用一对一路由（'c' 到 't'），但对于一对多、多对一或多对-很多。
# Select the node from which and to which the shortest path should be found.
from_node <- graph %>%
activate(nodes) %>%
filter(id_chr == "c") %>%
pull(id)
to_node <- graph %>%
activate(nodes) %>%
filter(id_chr == "t") %>%
pull(id)
# Find the shortest path between these nodes
path <- shortest_paths(
graph = graph,
from = from_node,
to = to_node,
output = 'both',
weights = graph %>% activate(edges) %>% pull(time)
)
生成的路径是一个列表，其中包含构成路径的节点和边。
$vpath
$vpath[[1]]
+ 18/20 vertices, from e43a089:
[1] 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
$epath
$epath[[1]]
+ 17/19 edges from e43a089:
[1] 3-- 4 4-- 5 5-- 6 6-- 7 7-- 8 8-- 9 9--10 10--11 11--12 12--13
[11] 13--14 14--15 15--16 16--17 17--18 18--19 19--20
我们可以创建一个只包含最短路径的节点和边的原始图的子图。
path_graph <- graph %>%
subgraph.edges(eids = path$epath %>% unlist()) %>%
as_tbl_graph()
# A tbl_graph: 18 nodes and 17 edges
#
# An undirected simple graph with 1 component
#
# Node Data: 18 x 4 (active)
id id_chr shape_id geometry
<int> <chr> <int> <POINT [°]>
1 3 c 52421 (-46.56212 -23.52124)
2 4 d 52421 (-46.56202 -23.52118)
3 5 e 52421 (-46.56193 -23.52113)
4 6 f 52421 (-46.56183 -23.52108)
5 7 g 52421 (-46.56181 -23.52121)
6 8 h 52421 (-46.56179 -23.52134)
# … with 12 more rows
#
# Edge Data: 17 x 6
from to length speed geometry time
<int> <int> [km] [km/h] <LINESTRING [°]> [h]
1 1 2 0.1914225 11 (-46.56212 -23.52124, -46.56202 -23.5… 0.017402…
2 2 3 0.1914225 11 (-46.56202 -23.52118, -46.56193 -23.5… 0.017402…
3 3 4 0.1914225 11 (-46.56193 -23.52113, -46.56183 -23.5… 0.017402…
# … with 14 more rows
在这里，发生了一些我不喜欢的事情。 Tidygraph/igraph 似乎有一个内部节点 ID 结构，你看到在子图中，egdes 数据中的 from 和 to 列与我们的 id 不匹配不再是节点数据中的列，而是简单地引用节点数据的行号。我不确定如何解决这个问题。
无论哪种方式，我们现在都有从 'c' 到 't' 的路径作为子图，并且可以轻松地对其进行分析。例如，通过计算路径的总行程时间（就像问题一样）。
path_graph %>%
activate(edges) %>%
as_tibble() %>%
summarise(total_time = sum(time))
# A tibble: 1 x 1
total_time
[h]
1 0.2958348
但绘制它也很容易，保留了地理信息（只需将节点和边导出为 sf 对象）。
ggplot() +
geom_sf(data = graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), col = 'darkgrey') +
geom_sf(data = graph %>% activate(nodes) %>% as_tibble() %>% st_as_sf(), col = 'darkgrey', size = 0.5) +
geom_sf(data = path_graph %>% activate(edges) %>% as_tibble() %>% st_as_sf(), lwd = 1, col = 'firebrick') +
geom_sf(data = path_graph %>% activate(nodes) %>% filter(id %in% c(from_node, to_node)) %>% as_tibble() %>% st_as_sf(), size = 2)
plot
关于这种 tidygraph-sf 方法的 r-spatial 博客 post 可能会出现 ;)