r - 如何将 sf 空间点列表转换为可路由图
问题描述
我有一个sf dataframe带有一系列点的对象，这些点代表公交路线的形状。我想把这个对象变成一个可路由的图，这样我就可以估计从点遍历c到t.
这是我尝试使用该dodgr 软件包的方法的方法，但我不确定我在这里做错了什么：
library(dodgr)
graph <- weight_streetnet(mydata, wt_profile = "motorcar", type_col="highway" , id_col = "id")
check_highway_osmid(x, wt_profile) 中的错误：请指定 type_col 用于加权 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"))
解决方案
如果您想将其包含在“整洁”的工作流程中，您还可以考虑使用 和 之间的sf混合tidygraph。后者为网络/图提供了一个整洁的框架，以tbl_graph类的形式，子类igraph（因此，您可以将tbl_graph所有igraph函数中的对象用作igraph对象）。但是，您可以将节点和边分析为 tibbles，并使用filter()、select()、mutate()等函数。当然，这些小标题也可以包含我们从中知道的几何列表列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
现在我们将所有内容都放在了一个图结构中，我们可以选择我们想要路由的节点，以及我们想要路由到的节点，并使用shortest_path函数 from以旅行时间作为权重变量找到它们之间的最短路径igraph。我们现在只使用一对一的路由（'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)
一篇关于这种 tidygraph-sf 方法的 r-spatial 博客文章可能会出现；）
推荐阅读
- python - 部署 django 项目后显示此错误 AttributeError: 'str' object has no attribute 'tzinfo'
- javascript - 从 D3.js 的 JSON 获取数据
- azure-storage - 将 $Logs 容器数据复制到另一个 Blob 位置
- parsing - 得到`PHP解析错误：语法错误，意外'：'，期待';' 或'{'`
- javascript - 验证以检查用户是否输入了正确的格式
- javascript - 多个 div 网格不响应阅读更多，阅读更少
- reactjs - Invariant Violation：试图注册两个同名的视图RNCSafeAreaProvider
- wordpress - 将域指向子目录
- flutter - 具有 .obs 和 hive 字段的颤振模型
- javascript - 如何从使用 Python 的网站中提取使用 javascript 生成的表数据？