r - 如何将 sf 空间点列表转换为可路由图

Question

我有一个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"))

score 2 · Accepted Answer

如果您想将其包含在“整洁”的工作流程中，您还可以考虑使用和之间的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 博客文章可能会出现；）

score 1 · Accepted Answer

我认为您可以通过将数据转换为igraph对象并使用igraph库中的功能来解决它。您需要建立边和顶点以及权重值。在igraph中，边是表示两个节点（源和目标）之间连接的链接。在这种情况下，链接是“街道”，点是节点。

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.

可以实现更复杂的网络和计算路线。例如，您可以设置道路的方向。

也许道奇可以解决这个问题，但我不确定。

score 1 · Accepted Answer

该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

请注意，对于像这样的“简单”问题，igraph通常会更快，因为它使用一组权重计算路线。在这种情况下，唯一真正的优势dodgr是能够使用“双权重” -$d_weighted和$d值 - 以便根据计算路线$d_weighted，并根据计算最终距离$d。

r - 如何将 sf 空间点列表转换为可路由图

可重现的数据

3 回答 3

Related

Reference