按路线类型提取 Openstreetmap 方式 ID 的交叉口计数
Extract count of intersections for Open Street Map way IDs, by route type
编辑并添加了更多详细信息
我有一个包含 2,061 个开放街道地图 (OSM) 路段的 shapefile。我的 shapefile 中的每个段都由其 OSM Way ID 标识。
这是我的数据中的五个细分的示例：
data =
structure(list(osm_id = structure(1:5, .Label = c("17990110",
"17993246", "17994983", "17994985", "17995338"), class = "factor"),
name = structure(c(1L, 3L, 4L, 5L, 2L), .Label = c("109th Avenue Northeast",
"85th Avenue Northeast", "Bunker Lake Boulevard Northeast",
"Foley Boulevard", "Northdale Boulevard Northwest"), class = "factor")), row.names = c(NA,
5L), class = c("sf", "data.frame"))
对于这 2061 个路段中的每一个，我想分别计算每种道路类型（住宅、主要、第三...）与其他道路的交叉点数量。
例如，this OSM way 与其他 27 条道路相交，其中 11 条是 "residential,"，3 条是 "secondary" 高速公路。
这对于分析来说是次要的，但最终，对于多种类型道路相交的交叉路口，我将 select "largest" 类型的道路。例如，this node连接一条服务道路和一条住宅道路；我想select小区路路。我想我可以为此创建自己的层次结构列表并在以后处理它。
我正在尝试使用 R Open Sci 包 osmdata
到目前为止，我可以使用 osmdata 到达 #2（信号交叉口）：
node_dat <- opq_osm_id(type = "node", id = '17990110')%>%
opq_string()%>%
osmdata_sf
node_dat_pts <- node_dat$osm_points
nrow(node_dat_pts[node_dat_pts$traffic_signals %in% "signal",])
这表明沿着该 OSM 路段有三个交通信号灯。 （尽管实际上只有两个信号交叉口——两个信号与分开的高速公路相关联……但这可能是另一个问题）。
我是 R 本地人，这就是 osmdata 包对我如此有吸引力的原因，但我愿意探索 Overpass API 中的查询。 TBH 我发现网站上的 example on how to get intersection nodes 不太适用——我不清楚如何将这个过程扩展到我拥有的 2000 多个方式 ID 的长列表（但如果文档或示例存在, 指给我看）。
我也愿意探索 Python 中的其他工具库； Python osmnx package 具有用于计算 "intersection density," 的优秀工具，但对于指定的多边形，如城市边界，似乎没有在 [=48= 内创建调用的功能] 通过方式或节点 ID。
我也知道我可以在 ArcGIS 或 QGIS 中做到这一点，但是因为我的数据库中已经有了这些 OSM ID，所以为一个大城市加载整个路口的该死 shapefile 似乎是一种浪费区域并进行某种缓冲过程以获取我需要的信息。另外，如果我有一个方便的脚本来通过方式或节点 ID 提取一些信息，我可以更轻松地扩大我提取的数据种类，以获得为 OSM 元素记录的重要信息的其他花絮。
感谢空间数据社区！
交通信号应该总是被标记"highway" = "traffic_signals"，但个别节点也可能被标记为"traffic_signals"。因此，获取所有交通信号的第一步可以这样完成：
library(osmdata)
hw <- opq("minneapolis") %>%
add_osm_feature(key = "highway") %>%
osmdata_sf()
signal_nodes <- opq("minneapolis") %>%
add_osm_feature(key = "traffic_signals") %>%
osmdata_sf()
index <- which (!is.na (hw$osm_points$traffic_signals) |
grepl ("traffic_signals", hw$osm_points$highway)) # grepl because can be "traffic_signals:<value>"
signal_node_ids <- unique (c (signal_nodes$osm_points$osm_id, hw$osm_points$osm_id [index]))
然后保存描述交通信号的节点的所有 OSM ID 值。
获得路口密度的一种直接方法是将高速公路的 sf 表示形式转换为 dodgr 网络，这是一个简单的 data.frame，每一行都是网络边缘。 poly2line 步骤将严格的 sf 多边形（例如环岛）转换为 linestring 对象，而 dodgr_contract_graph() 调用仅将网络减少到交汇点。
library(dodgr)
hw <- osm_poly2line(hw)$osm_lines %>%
weight_streetnet(keep_cols = "highway", wt_profile = 1) %>% # wt_profile irrelevant here
dodgr_contract_graph()
table(net$highway) 会给你不同种类高速公路的频率。然后，您可以按如下方式检查特定类型的结频率
net_footway <- net[net$highway == "footway", ]
table(table(c(net_footway$from_id, net_footway$to_id)))
值为1表示单向终端节点；值为 2 表示双向终端节点；值为 4 表示两条边之间的交叉点；等等。 table 上升到 14，因为人行道可能非常复杂，显然明尼阿波利斯某处有 7 条人行道的交汇处。这些 ID 是 OSM ID，因此您可以轻松检查它们中的哪些在 signal_node_ids 值中以确定哪些具有交通信号。
需要解决的剩余问题：
- 单独定义的
"highway" 类型之间的交集很容易，但不同类型之间的交集需要对此代码进行更复杂的修改。虽然简单明了，但您需要始终如一地对 dodgr data.frame 进行子集，其中边指向：$from_id -> $to_id.
- 将信号与特定路口相关联可能需要某种空间缓冲，因为正如您所建议的那样，单个路口可能有多个
"traffic_signals" 节点。请注意，没有 "right" 方法可以做到这一点，因为，例如，路口可能有单独的行人和汽车信号，是否将这些信号视为 "the same" 信号的决定总是需要一定程度的主观性。
我最终创建了一些依赖于 osmdat 包的自定义函数。我发现 osmdat 允许用户将自定义 API 调用传递给 Overpass。在 Overpass Turbo 中反复试验后，我想出了 Overpass 语法 "well enough" 来提取我需要的信息。我认为这三个独立的功能可以组合成一个 Overpass API 调用，但这超出了我的范围。
所以我首先创建了一个函数来获取与我的 "focal" 方式相关的所有方式的列表（我的数据框中的 2,061 个段中的一个）：
get_related_ways <- function(wayid, bboxstring){
related_ways_query <-
paste0("[bbox:", bboxstring,"];\n",
"way(id:", wayid, ");\n",
"rel(bw);\n", # get all sibling
"way(r);\n",
"out;")
related_ways <- osmdata_sf(related_ways_query)
related_ways <- related_ways$osm_lines
related_ways <- related_ways$osm_id
return(related_ways)
}
然后我创建了一个函数来为我的焦点方式提取节点 ID。
get_nodes_from_way <- function(wayid){
nodes_from_way <- opq_osm_id(id = wayid, "way")%>%osmdata_sf()
nodes_from_way <- nodes_from_way$osm_points
nodes_from_way$geometry <- NULL
setnames(nodes_from_way, old = 'osm_id', new = 'node_from_way_id')
return(nodes_from_way)
}
还有第三个函数，用于获取与我的焦点路径相交的所有路径的 ID。此函数需要输入焦点路径的节点 ID。
get_intersecting_ways <- function(nodeid, bboxstring){
node_to_intways_query <-
paste0("[bbox:", bboxstring,"];\n",
"node(id:", nodeid, ");\n",
"way(bn)[highway];\n",
"out;")
intways_from_node <- osmdata_sf(node_to_intways_query)
intways_from_node <- intways_from_node$osm_lines
intways_from_node$geometry <- NULL
intways_from_node <- intways_from_node[,names(intways_from_node) %in%c("osm_id", "name", "highway", "ref")]
setnames(intways_from_node, old = 'osm_id', new = 'intersecting_way_id')
setnames(intways_from_node, old = 'highway', new = 'intersecting_way_type')
setnames(intways_from_node, old = 'ref', new = 'intersecting_way_ref')
setnames(intways_from_node, old = 'name', new = 'intersecting_way_name')
return(intways_from_node)
}
对于这三个函数，我都有一个 "bboxstring" 或边界框字符串传递给 Overpass，希望能加快查询速度。哈。希望...
总之。我使用这三个函数创建了一个嵌套的 for 循环（不要评判我，我知道 purrr 存在，我只是觉得它们很直观！）。我还尝试通过使用 foreach 和 doParallel 将其并行化来砍掉我的方法，并将我的数据集分成 100 个块，每个块有 26 种方式。它仍然很慢。立交桥 API 可能很慢？更有可能是我在设置时做错了什么，这只是我第 3 次或第 4 次使用 doParallel。
for(this_part in unique(cmp_osmdat$partnum)){
osm_character_ids <- as.character(cmp_osmdat$osm_id)
# test:
# osm_character_ids <- osm_character_ids[1:3]
# for each parallel process to get our intersecting ways ("all ways")
all_ways <-
foreach(w = seq_along(osm_character_ids),
# require list of packages from above:
.packages = packs,
.errorhandling = "remove", # remove data frames that throw an error
# print the stuff:
.verbose = TRUE) %dopar% {
environmentIsLocked(asNamespace("curl"))
unlockBinding(sym = "has_internet", asNamespace("curl"))
assign(x = "has_internet", value = {function() T}, envir = asNamespace("curl"))
this_way_id <- osm_character_ids[[w]]
# find ways that are related to this one (same road, different segments)
# so that we can filter these out as intersections:
these_related_ways <- get_related_ways(this_way_id, this_bbox_string)
# get nodes of this way:
these_nodes_from_way <- get_nodes_from_way(this_way_id)
# adding a column to store this way id, for easy rbind later
# (foreach doesn't store list names?)
these_nodes_from_way$way_id <- this_way_id
# create an empty list to store interesecting ways for each node:
these_intersecting_ways <- list()
# get intersecing ways from nodes:
for(n in seq_along(these_nodes_from_way$node_from_way_id)){
this_node <- these_nodes_from_way$node_from_way_id[[n]]
# put intersecting ways into our empty list (the name of the list item will be the node ID)
these_intersecting_ways[[this_node]] <- get_intersecting_ways(this_node, this_bbox_string)
} # end get intersecting ways from node
# combine intersecting ways of each node into one data table:
these_intersecting_ways_df <- rbindlist(these_intersecting_ways, idcol = 'node_from_way_id', use.names = T, fill = T)
# get rid of intersections with this way's realtives (other segments of the same road):
these_intersecting_ways_df <- these_intersecting_ways_df[!these_intersecting_ways_df$intersecting_way_id %in% these_related_ways,]
# to get node information, merge intersecting ways to our node data:
nodes_and_ways <- merge(these_intersecting_ways_df, these_nodes_from_way, by = 'node_from_way_id')
# return node and intersection data
return(nodes_and_ways)
} # end foreach
nodes_and_ways_df <- rbindlist(all_ways, use.names = T, fill = T)
# save file, one for each part (results in 10 csvs)
write.csv(nodes_and_ways_df,
file = paste0("intersection_density_results/intersection-density-data-part-", this_part, ".csv"), row.names = F)
} # end 10 parts
stopCluster(cl)
这个过程的大体逻辑是：
- Select 来自数据帧块 (1-100) 的所有 WayID
- 从我的区块中的 26 种方式列表中获取一种方式 ID（"Focal way"）
- 找到焦点路径的路径 ID "siblings."
- 提取构成焦点路径的节点 ID（以及有关信号位置的信息 - TY、osmdata
- 对于焦点路径的每个节点id，找到与它相交的路径的路径id。也抢那些方式的分类。
- 摆脱实际上是 Focal Way 兄弟姐妹的 "intersecting ways" -- 这些部分是 Focal Way 的延续。 （例如，我会删除 this way from the list of intersecting ways if my Focal Way was this way
- 永远的 rbindlist
对于所有 2061 个片段，运行 大约需要 2-3 小时。那是很长的时间；但是，即使在 Overpass Turbo 中直接查询也很慢，所以...也许这是正确的。
编辑并添加了更多详细信息
我有一个包含 2,061 个开放街道地图 (OSM) 路段的 shapefile。我的 shapefile 中的每个段都由其 OSM Way ID 标识。
这是我的数据中的五个细分的示例：
data =
structure(list(osm_id = structure(1:5, .Label = c("17990110",
"17993246", "17994983", "17994985", "17995338"), class = "factor"),
name = structure(c(1L, 3L, 4L, 5L, 2L), .Label = c("109th Avenue Northeast",
"85th Avenue Northeast", "Bunker Lake Boulevard Northeast",
"Foley Boulevard", "Northdale Boulevard Northwest"), class = "factor")), row.names = c(NA,
5L), class = c("sf", "data.frame"))
对于这 2061 个路段中的每一个，我想分别计算每种道路类型（住宅、主要、第三...）与其他道路的交叉点数量。
例如，this OSM way 与其他 27 条道路相交，其中 11 条是 "residential,"，3 条是 "secondary" 高速公路。
这对于分析来说是次要的，但最终，对于多种类型道路相交的交叉路口，我将 select "largest" 类型的道路。例如，this node连接一条服务道路和一条住宅道路；我想select小区路路。我想我可以为此创建自己的层次结构列表并在以后处理它。
我正在尝试使用 R Open Sci 包 osmdata
到目前为止，我可以使用 osmdata 到达 #2（信号交叉口）：
node_dat <- opq_osm_id(type = "node", id = '17990110')%>%
opq_string()%>%
osmdata_sf
node_dat_pts <- node_dat$osm_points
nrow(node_dat_pts[node_dat_pts$traffic_signals %in% "signal",])
这表明沿着该 OSM 路段有三个交通信号灯。 （尽管实际上只有两个信号交叉口——两个信号与分开的高速公路相关联……但这可能是另一个问题）。
我是 R 本地人，这就是 osmdata 包对我如此有吸引力的原因，但我愿意探索 Overpass API 中的查询。 TBH 我发现网站上的 example on how to get intersection nodes 不太适用——我不清楚如何将这个过程扩展到我拥有的 2000 多个方式 ID 的长列表（但如果文档或示例存在, 指给我看）。
我也愿意探索 Python 中的其他工具库； Python osmnx package 具有用于计算 "intersection density," 的优秀工具，但对于指定的多边形，如城市边界，似乎没有在 [=48= 内创建调用的功能] 通过方式或节点 ID。
我也知道我可以在 ArcGIS 或 QGIS 中做到这一点，但是因为我的数据库中已经有了这些 OSM ID，所以为一个大城市加载整个路口的该死 shapefile 似乎是一种浪费区域并进行某种缓冲过程以获取我需要的信息。另外，如果我有一个方便的脚本来通过方式或节点 ID 提取一些信息，我可以更轻松地扩大我提取的数据种类，以获得为 OSM 元素记录的重要信息的其他花絮。
感谢空间数据社区！
交通信号应该总是被标记"highway" = "traffic_signals"，但个别节点也可能被标记为"traffic_signals"。因此，获取所有交通信号的第一步可以这样完成：
library(osmdata)
hw <- opq("minneapolis") %>%
add_osm_feature(key = "highway") %>%
osmdata_sf()
signal_nodes <- opq("minneapolis") %>%
add_osm_feature(key = "traffic_signals") %>%
osmdata_sf()
index <- which (!is.na (hw$osm_points$traffic_signals) |
grepl ("traffic_signals", hw$osm_points$highway)) # grepl because can be "traffic_signals:<value>"
signal_node_ids <- unique (c (signal_nodes$osm_points$osm_id, hw$osm_points$osm_id [index]))
然后保存描述交通信号的节点的所有 OSM ID 值。
获得路口密度的一种直接方法是将高速公路的 sf 表示形式转换为 dodgr 网络，这是一个简单的 data.frame，每一行都是网络边缘。 poly2line 步骤将严格的 sf 多边形（例如环岛）转换为 linestring 对象，而 dodgr_contract_graph() 调用仅将网络减少到交汇点。
library(dodgr)
hw <- osm_poly2line(hw)$osm_lines %>%
weight_streetnet(keep_cols = "highway", wt_profile = 1) %>% # wt_profile irrelevant here
dodgr_contract_graph()
table(net$highway) 会给你不同种类高速公路的频率。然后，您可以按如下方式检查特定类型的结频率
net_footway <- net[net$highway == "footway", ]
table(table(c(net_footway$from_id, net_footway$to_id)))
值为1表示单向终端节点；值为 2 表示双向终端节点；值为 4 表示两条边之间的交叉点；等等。 table 上升到 14，因为人行道可能非常复杂，显然明尼阿波利斯某处有 7 条人行道的交汇处。这些 ID 是 OSM ID，因此您可以轻松检查它们中的哪些在 signal_node_ids 值中以确定哪些具有交通信号。
需要解决的剩余问题：
- 单独定义的
"highway"类型之间的交集很容易，但不同类型之间的交集需要对此代码进行更复杂的修改。虽然简单明了，但您需要始终如一地对dodgr data.frame进行子集，其中边指向：$from_id -> $to_id. - 将信号与特定路口相关联可能需要某种空间缓冲，因为正如您所建议的那样，单个路口可能有多个
"traffic_signals"节点。请注意，没有 "right" 方法可以做到这一点，因为，例如，路口可能有单独的行人和汽车信号，是否将这些信号视为 "the same" 信号的决定总是需要一定程度的主观性。
我最终创建了一些依赖于 osmdat 包的自定义函数。我发现 osmdat 允许用户将自定义 API 调用传递给 Overpass。在 Overpass Turbo 中反复试验后，我想出了 Overpass 语法 "well enough" 来提取我需要的信息。我认为这三个独立的功能可以组合成一个 Overpass API 调用，但这超出了我的范围。
所以我首先创建了一个函数来获取与我的 "focal" 方式相关的所有方式的列表（我的数据框中的 2,061 个段中的一个）：
get_related_ways <- function(wayid, bboxstring){
related_ways_query <-
paste0("[bbox:", bboxstring,"];\n",
"way(id:", wayid, ");\n",
"rel(bw);\n", # get all sibling
"way(r);\n",
"out;")
related_ways <- osmdata_sf(related_ways_query)
related_ways <- related_ways$osm_lines
related_ways <- related_ways$osm_id
return(related_ways)
}
然后我创建了一个函数来为我的焦点方式提取节点 ID。
get_nodes_from_way <- function(wayid){
nodes_from_way <- opq_osm_id(id = wayid, "way")%>%osmdata_sf()
nodes_from_way <- nodes_from_way$osm_points
nodes_from_way$geometry <- NULL
setnames(nodes_from_way, old = 'osm_id', new = 'node_from_way_id')
return(nodes_from_way)
}
还有第三个函数，用于获取与我的焦点路径相交的所有路径的 ID。此函数需要输入焦点路径的节点 ID。
get_intersecting_ways <- function(nodeid, bboxstring){
node_to_intways_query <-
paste0("[bbox:", bboxstring,"];\n",
"node(id:", nodeid, ");\n",
"way(bn)[highway];\n",
"out;")
intways_from_node <- osmdata_sf(node_to_intways_query)
intways_from_node <- intways_from_node$osm_lines
intways_from_node$geometry <- NULL
intways_from_node <- intways_from_node[,names(intways_from_node) %in%c("osm_id", "name", "highway", "ref")]
setnames(intways_from_node, old = 'osm_id', new = 'intersecting_way_id')
setnames(intways_from_node, old = 'highway', new = 'intersecting_way_type')
setnames(intways_from_node, old = 'ref', new = 'intersecting_way_ref')
setnames(intways_from_node, old = 'name', new = 'intersecting_way_name')
return(intways_from_node)
}
对于这三个函数，我都有一个 "bboxstring" 或边界框字符串传递给 Overpass，希望能加快查询速度。哈。希望...
总之。我使用这三个函数创建了一个嵌套的 for 循环（不要评判我，我知道 purrr 存在，我只是觉得它们很直观！）。我还尝试通过使用 foreach 和 doParallel 将其并行化来砍掉我的方法，并将我的数据集分成 100 个块，每个块有 26 种方式。它仍然很慢。立交桥 API 可能很慢？更有可能是我在设置时做错了什么，这只是我第 3 次或第 4 次使用 doParallel。
for(this_part in unique(cmp_osmdat$partnum)){
osm_character_ids <- as.character(cmp_osmdat$osm_id)
# test:
# osm_character_ids <- osm_character_ids[1:3]
# for each parallel process to get our intersecting ways ("all ways")
all_ways <-
foreach(w = seq_along(osm_character_ids),
# require list of packages from above:
.packages = packs,
.errorhandling = "remove", # remove data frames that throw an error
# print the stuff:
.verbose = TRUE) %dopar% {
environmentIsLocked(asNamespace("curl"))
unlockBinding(sym = "has_internet", asNamespace("curl"))
assign(x = "has_internet", value = {function() T}, envir = asNamespace("curl"))
this_way_id <- osm_character_ids[[w]]
# find ways that are related to this one (same road, different segments)
# so that we can filter these out as intersections:
these_related_ways <- get_related_ways(this_way_id, this_bbox_string)
# get nodes of this way:
these_nodes_from_way <- get_nodes_from_way(this_way_id)
# adding a column to store this way id, for easy rbind later
# (foreach doesn't store list names?)
these_nodes_from_way$way_id <- this_way_id
# create an empty list to store interesecting ways for each node:
these_intersecting_ways <- list()
# get intersecing ways from nodes:
for(n in seq_along(these_nodes_from_way$node_from_way_id)){
this_node <- these_nodes_from_way$node_from_way_id[[n]]
# put intersecting ways into our empty list (the name of the list item will be the node ID)
these_intersecting_ways[[this_node]] <- get_intersecting_ways(this_node, this_bbox_string)
} # end get intersecting ways from node
# combine intersecting ways of each node into one data table:
these_intersecting_ways_df <- rbindlist(these_intersecting_ways, idcol = 'node_from_way_id', use.names = T, fill = T)
# get rid of intersections with this way's realtives (other segments of the same road):
these_intersecting_ways_df <- these_intersecting_ways_df[!these_intersecting_ways_df$intersecting_way_id %in% these_related_ways,]
# to get node information, merge intersecting ways to our node data:
nodes_and_ways <- merge(these_intersecting_ways_df, these_nodes_from_way, by = 'node_from_way_id')
# return node and intersection data
return(nodes_and_ways)
} # end foreach
nodes_and_ways_df <- rbindlist(all_ways, use.names = T, fill = T)
# save file, one for each part (results in 10 csvs)
write.csv(nodes_and_ways_df,
file = paste0("intersection_density_results/intersection-density-data-part-", this_part, ".csv"), row.names = F)
} # end 10 parts
stopCluster(cl)
这个过程的大体逻辑是：
- Select 来自数据帧块 (1-100) 的所有 WayID
- 从我的区块中的 26 种方式列表中获取一种方式 ID（"Focal way"）
- 找到焦点路径的路径 ID "siblings."
- 提取构成焦点路径的节点 ID（以及有关信号位置的信息 - TY、osmdata
- 对于焦点路径的每个节点id，找到与它相交的路径的路径id。也抢那些方式的分类。
- 摆脱实际上是 Focal Way 兄弟姐妹的 "intersecting ways" -- 这些部分是 Focal Way 的延续。 （例如，我会删除 this way from the list of intersecting ways if my Focal Way was this way
- 永远的 rbindlist
对于所有 2061 个片段，运行 大约需要 2-3 小时。那是很长的时间；但是，即使在 Overpass Turbo 中直接查询也很慢，所以...也许这是正确的。