R:如何将 osmdata 的运行时间减少到 igraph 转换
R: How could I reduce runtime of osmdata to an igraph conversion
是否可以减少以下代码的运行时间?
我的目标是从由框边界指定的开放街道数据区域获取加权 igraph 对象。
目前我正在尝试使用立交桥 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 对象到具有函数 [=15= 的 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)
中缺少此内容
如果这可以添加到 silicate 到 igraph 转换
我可以跳过 dat_osmar <-as_osmar(xmlParse(dat))
步骤
然后走 overpass->silicate->igraph
比 overpass->xml->osmar->igraph
.
快得多的路线
osmdata 包还提供 sf 响应 osmdata_sf()
所以也许工作流程 overpass->sf->igraph
更快,但在使用这种方式时,我需要根据边缘距离将权重合并到图中,我目前还不够好,我会非常感激任何帮助。
此外,在使用 sf 和生成的 igraph 对象时,不应丢失 openstreetmap gps 点及其 ID 之间的连接。这意味着我应该能够从生成的 Igraph 中找到一个 ID 的 gps 位置。查找 table 就足够了。如果我走 overpass->silicate->igraph
或 overpass->xml->osmar->igraph
路线,这是可能的。我不确定 overpass->sf->igraph
路线是否仍然可行。
似乎将 xml 数据转换为另一种格式需要很长时间。与其使用 xml,不如向 return 请求立交桥 sf
对象并使用它可能会更快。然后 sf
对象可以被 sfnetworks
包操作和使用以获得网络,同时保留网络的空间方面。 sfnetworks
(或tidygraph
)包中的函数可以添加权重。
我认为下面的代码着重于处理速度和边缘权重问题。其他问题,如寻找最近的节点或边缘,可以使用 sf
包的功能来解决,但没有解决。否则这不仅仅是一个一次性的 SO 问题。
通过对边缘使用 st_simplify
,可以提高速度,但会牺牲空间精度。这种方法的一个问题是 stnetworks 在每个线串与另一个线串相遇的地方放置了一个节点。数据 returned 通常将一条道路分成多条线串。例如,在下面的边缘图中可以看到两条较长的黄色道路。可能是一个可以解决的问题,但在这种情况下可能不值得。
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
只有边的图和有节点的边的图:
几条长路歪斜了颜色,但说明了一条路一分为二。
编辑:
通过纬度/经度坐标将注释添加到 select 最近的边缘(道路)。节点(上面的交叉点/红点)没有我知道的 osm id 号。节点由 sfnetworks
.
创建
从 lat/lon 点的 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
的作者,也许他会有一些见解。
如果您想在 R 中创建一个从道路网络开始的图形对象,那么我将使用以下过程。
首先,我需要从 github 存储库安装 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 and 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-03-26 创建
如果您想阅读有关 sfnetworks
的更多详细信息,则可以查看 website and the introductory vignettes。话虽这么说,我不明白你的意思
connection between openstreetmap gps points and their IDs should not be lost
您能否通过对原始问题的评论或编辑来添加更多详细信息?为什么需要 OSM id? OSM id 是什么意思?我想我需要更多细节来扩展这个答案。
编辑
我刚刚重读了@mrhellmann 的回答,我注意到我忘记了将 POLYGON 数据转换为线。无论如何,我建议在 运行 之后立即应用 osmdata::osm_poly2line()
通过立交桥 API.
下载 OSM 数据的代码
是否可以减少以下代码的运行时间?
我的目标是从由框边界指定的开放街道数据区域获取加权 igraph 对象。
目前我正在尝试使用立交桥 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 对象到具有函数 [=15= 的 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)
如果这可以添加到 silicate 到 igraph 转换
我可以跳过 dat_osmar <-as_osmar(xmlParse(dat))
步骤
然后走 overpass->silicate->igraph
比 overpass->xml->osmar->igraph
.
osmdata 包还提供 sf 响应 osmdata_sf()
所以也许工作流程 overpass->sf->igraph
更快,但在使用这种方式时,我需要根据边缘距离将权重合并到图中,我目前还不够好,我会非常感激任何帮助。
此外,在使用 sf 和生成的 igraph 对象时,不应丢失 openstreetmap gps 点及其 ID 之间的连接。这意味着我应该能够从生成的 Igraph 中找到一个 ID 的 gps 位置。查找 table 就足够了。如果我走 overpass->silicate->igraph
或 overpass->xml->osmar->igraph
路线,这是可能的。我不确定 overpass->sf->igraph
路线是否仍然可行。
似乎将 xml 数据转换为另一种格式需要很长时间。与其使用 xml,不如向 return 请求立交桥 sf
对象并使用它可能会更快。然后 sf
对象可以被 sfnetworks
包操作和使用以获得网络,同时保留网络的空间方面。 sfnetworks
(或tidygraph
)包中的函数可以添加权重。
我认为下面的代码着重于处理速度和边缘权重问题。其他问题,如寻找最近的节点或边缘,可以使用 sf
包的功能来解决,但没有解决。否则这不仅仅是一个一次性的 SO 问题。
通过对边缘使用 st_simplify
,可以提高速度,但会牺牲空间精度。这种方法的一个问题是 stnetworks 在每个线串与另一个线串相遇的地方放置了一个节点。数据 returned 通常将一条道路分成多条线串。例如,在下面的边缘图中可以看到两条较长的黄色道路。可能是一个可以解决的问题,但在这种情况下可能不值得。
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
只有边的图和有节点的边的图:
几条长路歪斜了颜色,但说明了一条路一分为二。
编辑:
通过纬度/经度坐标将注释添加到 select 最近的边缘(道路)。节点(上面的交叉点/红点)没有我知道的 osm id 号。节点由 sfnetworks
.
从 lat/lon 点的 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
的作者,也许他会有一些见解。
如果您想在 R 中创建一个从道路网络开始的图形对象,那么我将使用以下过程。
首先,我需要从 github 存储库安装 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 and 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-03-26 创建如果您想阅读有关 sfnetworks
的更多详细信息,则可以查看 website and the introductory vignettes。话虽这么说,我不明白你的意思
connection between openstreetmap gps points and their IDs should not be lost
您能否通过对原始问题的评论或编辑来添加更多详细信息?为什么需要 OSM id? OSM id 是什么意思?我想我需要更多细节来扩展这个答案。
编辑
我刚刚重读了@mrhellmann 的回答,我注意到我忘记了将 POLYGON 数据转换为线。无论如何,我建议在 运行 之后立即应用 osmdata::osm_poly2line()
通过立交桥 API.