R：如何将 osmdata 的运行时间缩短为 igraph 转换答案

【问题标题】：R: How could I reduce runtime of osmdata to an igraph conversionR：如何将 osmdata 的运行时间缩短为 igraph 转换
【发布时间】：2021-06-19 23:13:36
【问题描述】：

是否可以减少以下代码的运行时间？

我的目标是从框边界指定的开放街道数据区域中获取加权 igraph 对象。

目前我正在尝试使用 overpass api 来减轻内存负载，因此我不必在内存中保留大的 osm 文件。

首先我得到一个由 bbox（仅街道）指定的 osm 数据作为 xml 结构

library(osmdata)
library(osmar)
install.packages("remotes")
remotes::install_github("hypertidy/scgraph")
library(scgraph)


dat <- opq(bbox = c(11.68771, 47.75233, 12.35058, 48.19743 )) %>%
  add_osm_feature(key = 'highway',value = c("trunk", "trunk_link", "primary","primary_link", "secondary", "secondary_link", "tertiary","tertiary_link", "residential", "unclassified" ))%>% 
  osmdata_xml ()

然后我将生成的 xml 对象 dat 转换为 osmar 对象 dat_osmar，最后转换为 igraph 对象：

dat_osmar <-as_osmar(xmlParse(dat))
dat_graoh <- as_igraph(dat_osmar)

如何优化这些例程？

也许可以将 dat (XML) 对象分成块并并行解析？

我经过几个步骤才最终得到一个加权无向图。

目前整个过程在我的机器上需要 89.555 秒。

如果我可以缩短这两个步骤的运行时间：

dat_osmar <-as_osmar(xmlParse(dat))
dat_graoh <- as_igraph(dat_osmar)

这已经有帮助了。

我尝试的一种方法是使用 osmdata_sc() 而不是 osmdata_xml()。

这提供了一个硅酸盐对象，我可以将其转换为：

scgraph::sc_as_igraph(dat)

到 igraph。

它相当快，但遗憾的是重量正在丢失，所以它不是一个解决方案。

原因是：如果我使用从 osmar 对象到具有函数osmar::as_igraph() 的 igraph 对象的转换，则权重是根据两者之间的距离计算的两条边并添加到 igraph：

    edges <- lapply(dat, function(x) {
    n <- nrow(x)
    from <- 1:(n - 1)
    to <- 2:n
    weights <- distHaversine(x[from, c("lon", "lat")], x[to, 
      c("lon", "lat")])
    cbind(from_node_id = x[from, "ref"], to_node_id = x[to, 
      "ref"], way_id = x[1, "id"], weights = weights)
  })

scgraph::sc_as_igraph(dat) 中缺少此内容

如果这可以添加到 硅酸盐 到 igraph 转换我可以跳过dat_osmar <-as_osmar(xmlParse(dat)) 步骤然后走overpass->silicate->igraph 路线，它比overpass->xml->osmar->igraph 快得多。

osmdata 包还通过 osmdata_sf()

提供 sf 响应

所以也许overpass->sf->igraph 的工作流程更快，但在使用这种方式时，我需要根据边的距离将权重合并到图表中，而我目前还不够好，非常感谢任何帮助。

此外，在使用 sf 和生成的 igraph 对象时，openstreetmap gps 点与其 ID 之间的连接不应丢失。这意味着我应该能够从生成的 Igraph 中找到一个 ID 的 gps 位置。一个查找表就足够了。如果我去overpass->silicate->igraph 或overpass->xml->osmar->igraph 路线，这是可能的。我不确定overpass->sf->igraph 路由是否仍然可行。

【问题讨论】：

嗨！如果您对osmar 的替代方法感兴趣，我可以尝试提供一个基于名为sfnetworks 的R 包的解决方案。 sfnetworks是基于tidygraph，这意味着sfnetwork返回的对象也是igraph对象。
您好，是的。我需要能够做两件事：1. 根据 bbox 区域从立交桥获取数据，2. 从该区域获取 igraph 3. 能够将 igraph 连接/查找到 openstreetmap ID 到 gps 位置。例如，现在我使用 igraph 进行路由/图形处理，并在 osmar 对象中查找 ID 和 gps 值以及 openstreetmap 数据的其他键。如果可以通过 sfnetworks 和快速实现相同的效果（如果不超过 10 秒，不包括访问立交桥会很棒），我会非常乐意使用这个解决方案！

标签： r openstreetmap igraph sf overpass-api

【解决方案1】：

如果您想从 R 中的道路网络创建一个图形对象，那么我将使用以下过程。

首先，我需要从 github repo 安装sfnetworks（因为我们最近修复了一些错误并且最新版本不在 CRAN 上）

remotes::install_github("luukvdmeer/sfnetworks", quiet = TRUE)

然后加载包

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(sfnetworks)
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

从 Overpass API 下载数据

my_osm_data <- opq(bbox = c(11.68771, 47.75233, 12.35058, 48.19743 )) %>%
  add_osm_feature(
    key = 'highway', 
    value = c("trunk", "trunk_link", "primary","primary_link", "secondary", "secondary_link", "tertiary","tertiary_link", "residential", "unclassified")
  ) %>% 
  osmdata_sf(quiet = FALSE)
#> Issuing query to Overpass API ...
#> Rate limit: 2
#> Query complete!
#> converting OSM data to sf format

现在我提取道路并构建 sfnetwork 对象：

system.time({
  # extract the roads
  my_roads <- st_geometry(my_osm_data$osm_lines)
  
  # build the sfnetwork object
  my_sfn <- as_sfnetwork(my_roads, directed = FALSE, length_as_weight = TRUE)
})
#>    user  system elapsed 
#>    3.03    0.16    3.28

如您所见，下载 OSM 数据后，运行该过程只需几秒钟。

目前我忽略了my_osm_data$osm_lines中的所有字段，但是如果您需要将my_osm_data$osm_lines中的某些列添加到my_roads中，那么您可以将之前的代码修改为：my_roads <- my_osm_data$osm_lines[, "relevant columns"]。

关于构造sfnetwork 对象的一些细节：参数"directed = FALSE" 指定我们要构建一个无向图（有关详细信息，请参阅文档here 和here），而参数length_as_weight = TRUE 表示边的长度将存储在名为 "weight" 的列中，并由 igraph 和 tidygraph 算法使用。

这是my_sfn对象的打印：

my_sfn
#> # A sfnetwork with 33179 nodes and 28439 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # An undirected multigraph with 6312 components with spatially explicit edges
#> #
#> Registered S3 method overwritten by 'cli':
#>   method     from         
#>   print.boxx spatstat.geom
#> # Node Data:     33,179 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 11.6757 ymin: 47.74745 xmax: 12.39161 ymax: 48.22025
#>                     x
#>           <POINT [°]>
#> 1 (11.68861 47.90971)
#> 2 (11.68454 47.90937)
#> 3 (11.75216 48.17638)
#> 4 (11.75358 48.17438)
#> 5  (11.7528 48.17351)
#> 6 (11.74822 48.17286)
#> # ... with 33,173 more rows
#> #
#> # Edge Data:     28,439 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 11.6757 ymin: 47.74745 xmax: 12.39161 ymax: 48.22025
#>    from    to                                                           x weight
#>   <int> <int>                                            <LINESTRING [°]>  <dbl>
#> 1     1     2 (11.68861 47.90971, 11.6878 47.90965, 11.68653 47.90954, 1~   306.
#> 2     3     4 (11.75216 48.17638, 11.75224 48.17626, 11.75272 48.17556, ~   246.
#> 3     5     6 (11.7528 48.17351, 11.75264 48.17344, 11.75227 48.17329, 1~   382.
#> # ... with 28,436 more rows

my_sfn 根据定义是一个 igraph 对象：

class(my_sfn)
#> [1] "sfnetwork" "tbl_graph" "igraph"

但是，如果你想更明确，那么

as.igraph(my_sfn)
#> IGRAPH 101dcdf U-W- 33179 28439 -- 
#> + attr: x (v/x), x (e/x), weight (e/n)
#> + edges from 101dcdf:
#>  [1]   1--  2   3--  4   5--  6   7--  8   9-- 10  11-- 12  13-- 14  15-- 16
#>  [9]  17-- 18  16-- 19  20-- 21  21-- 22  23-- 24  25-- 26  27-- 28  29-- 30
#> [17]  31-- 32  33-- 34  35-- 36  37-- 38  39-- 40  41-- 42  43-- 44  45-- 46
#> [25]  14-- 47  48-- 49  50-- 51  52-- 53  54-- 55  56-- 57  36-- 58  58-- 59
#> [33]  60-- 61  62-- 63  64-- 65  66-- 67  68-- 69  70-- 71  72-- 73  74-- 75
#> [41]  76-- 77  78-- 79  80-- 81  82-- 83  84-- 85  86-- 87  88-- 89  90-- 91
#> [49]  92-- 93  94-- 95  96-- 97  98-- 99 100--101 102--103 104--105 106--107
#> [57] 108--109 110--111 112--113  80--114 115--116 117--118 119--120 121--122
#> + ... omitted several edges

可以看到边有一个权重属性，等于每个 LINESTRING 几何体的长度：

all.equal(
  target = igraph::edge_attr(as.igraph(my_sfn), "weight"), 
  current = as.numeric(st_length(my_roads))
)
#> [1] TRUE

^{由reprex package (v1.0.0) 于 2021 年 3 月 26 日创建}

如果您想阅读有关sfnetworks 的更多详细信息，则可以查看website 和introductory vignettes。话虽如此，我不明白你的意思

openstreetmap gps 点与其 ID 之间的连接不应丢失

您能否通过评论或对原始问题的编辑添加更多详细信息？为什么需要 OSM id？ OSM id 是什么意思？我认为我需要更多细节来扩展这个答案。

编辑

我刚刚重新阅读了@mrhellmann 的答案，我注意到我忘记将 POLYGON 数据转换为线。无论如何，我建议在运行代码后立即申请 osmdata::osm_poly2line() 通过 Overpass API 下载 OSM 数据。

【讨论】：

【解决方案2】：

似乎将 xml 数据转换为另一种格式需要很长时间。而不是使用 xml，而是要求 overpass 返回一个 sf 对象并使用它可能会更快。然后sf 对象可以被sfnetworks 包操作和使用以获取网络，同时保留网络的空间方面。权重可以通过sfnetworks（或tidygraph）包中的函数添加。

我认为下面的代码侧重于处理速度和边缘权重问题。其他问题，例如查找最近的节点或边，可以使用sf 包的功能解决，但未解决。否则，这不仅仅是一个一次性的 SO 问题。

通过将st_simplify 用于边缘，可以以空间精度为代价提高速度。这种方法的一个问题是 stnetworks 会在每个线串与另一个线串相交的地方放置一个节点。返回的数据通常将一条道路分割成多个线串。例如，在下面的边图中看到两条较长的黄色道路。可能是一个可以解决的问题，但在这种情况下可能不值得。

library(osmdata)
#library(osmar)
library(tidyverse)
library(sf)
library(ggplot2)
library(sfnetworks)
library(tidygraph)


# get data as an sf object rather than xml
## This is the slowest part of the code.
dat_sf <- opq(bbox = c(11.68771, 47.75233, 12.35058, 48.19743 )) %>%
  add_osm_feature(key = 'highway',value = c("trunk", "trunk_link", "primary","primary_link", "secondary", "secondary_link", "tertiary","tertiary_link", "residential", "unclassified" ))%>% 
  osmdata_sf()



# Only keep lines & polygons. Points take up too much memory & 
## all seem to be on lines anyway. Change polygons to LINESTRING,
## as most seem to be roundabouts or a few odd cases.
lines_sf <- dat_sf$osm_lines %>% select(osm_id) %>% st_sf()
polys_sf <- dat_sf$osm_polygons %>% select(osm_id) %>% st_sf() %>% 
  st_cast('LINESTRING')

# Combine the two above sf objects into one 
dat_sf_bound <- rbind(lines_sf, polys_sf)

# get an sfnetwork
dat_sf_net <- as_sfnetwork(dat_sf_bound)

# add edge weight as distance
dat_sf_net <- dat_sf_net %>%
  activate(edges) %>%
  mutate(weight = edge_length())

dat_sf_net 应该是这样的：

> dat_sf_net
# An sfnetwork with 33255 nodes and 28608 edges
#
# CRS:  EPSG:4326 
#
# A directed multigraph with 6391 components with spatially explicit edges
#
# Edge Data:     28,608 x 4 (active)
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: 11.6757 ymin: 47.74745 xmax: 12.39161 ymax: 48.22025
   from    to   weight                                                                  geometry
  <int> <int>      [m]                                                          <LINESTRING [°]>
1     1     2 306.3998 (11.68861 47.90971, 11.6878 47.90965, 11.68653 47.90954, 11.68597 47.909…
2     3     4 245.9225 (11.75216 48.17638, 11.75224 48.17626, 11.75272 48.17556, 11.7528 48.175…
3     5     6 382.2423 (11.7528 48.17351, 11.75264 48.17344, 11.75227 48.17329, 11.752 48.1732,…
4     7     8 131.1373 (11.70029 47.87861, 11.70046 47.87869, 11.70069 47.87879, 11.70138 47.87…
5     9    10 252.9170 (11.75733 48.17339, 11.75732 48.17343, 11.75726 48.17357, 11.75718 48.17…
6    11    12 131.6942 (11.75582 48.17036, 11.75551 48.1707, 11.75521 48.17106, 11.75507 48.171…
# … with 28,602 more rows
#
# Node Data:     33,255 x 1
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 11.6757 ymin: 47.74745 xmax: 12.39161 ymax: 48.22025
             geometry
          <POINT [°]>
1 (11.68861 47.90971)
2 (11.68454 47.90937)
3 (11.75216 48.17638)
# … with 33,252 more rows

仅边的图，以及带有节点的边的图：

几条长长的道路歪斜了颜色，但说明了一条分成两半的道路。

编辑：解决注释以通过纬度/经度坐标选择最近的边缘（道路）。节点（上面的交叉点/红点）没有我知道的 osm id 编号。节点由sfnetworks创建。

从一个sf 对象的纬度/经度点开始作为我们组成的gps 坐标。

# random point
gps <- sfheaders::sf_point(data.frame(x = 11.81854, y = 48.04514)) %>% st_set_crs(4326)

# nearest edge(road) to the point. dat_sf_net must have edges activated.
near_edge <- st_nearest_feature(gps, dat_sf_net %>% st_as_sf())

>near_edge
[1] 4359

> st_as_sf(dat_sf_net)[near_edge,]
Simple feature collection with 1 feature and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 11.81119 ymin: 48.02841 xmax: 11.82061 ymax: 48.06845
Geodetic CRS:  WGS 84
# A tibble: 1 x 5
   from    to osm_id                                                           geometry   weight
  <int> <int> <chr>                                                    <LINESTRING [°]>      [m]
1  7590  7591 24232418 (11.81289 48.02841, 11.81213 48.03014, 11.81202 48.03062, 11.81… 4576.273

p3 <- gplot() +
      geom_sf(data = st_as_sf(dat_sf_net), color = 'black') +       
      geom_sf(data = gps, color = 'red') + 
      geom_sf(data = st_as_sf(dat_sf_net)[near_edge,], color = 'orange') + 
      coord_sf(xlim = c(11.7, 11.9), ylim = c(48, 48.1))

看起来像上面评论的@agila。既然是sfnetworks的作者，说不定会有一些见解。

【讨论】：

谢谢！我将如何查询数据以获取 gps 位置？目前，我使用 gps 位置查询 osm 数据对象，以在 osm 数据中找到一个近似点，以获取 igraph 中正在使用的靠近它的节点 ID。然后我拿 id 并且可以在 igraph 中找到这个 id。由于我不再有 osm 对象，我将如何在此解决方案中进行处理？我可以直接查询 dat_sf_net 以获取近似 gps 点还是 sf 对象是否包含到 openstreetmap 节点 ID 连接的 gps？
步骤：near_edge % st_as_sf()) 抛出错误：错误在st_nearest_feature(gps, dat_sf_net %>% st_as_sf()) : st_crs(x) = = st_crs(y) 不正确