Edited with additional details added
I have a shapefile of 2,061 Open Street Map (OSM) road segments. Each segment in my shapefile is is identified by its OSM Way ID.
Here is an example of five segments from my data:
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"))
For each of these 2061 road segments, I want to count the number of intersections with other roads, separately for each road type (residential, primary, tertiary...).
For example, this OSM way intersects 27 other ways, 11 of which are "residential," and 3 of which are "secondary" highways.
This is secondary to the analysis, but ultimately, for intersections where multiple types of roads join, I will select the "largest" type of road. For example, this node joins a service road and residential road; I would like to select the residential road road. I think I can create my own hierarchy list for this and deal with it later.
I am trying to use R Open Sci package osmdata
So far I can get at #2 (signalized intersections) using osmdata:
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",])
This reveals that there are three traffic signals along this OSM segment. (Although, in reality, there are only two signalized intersections -- two signals are associated with a divided highway...but that might be a problem for another time).
I'm an R native, which is why the osmdata package is so appealing to me, but I'm open to exploring querying in Overpass API. TBH I found the example on how to get intersection nodes on the website not quite applicable -- and I'm unclear of how to scale up this process to the long list of 2000+ way IDs that I have (but if the docs or an example exist, point me to it).
I'm also open to exploring other tool libraries in Python; the Python osmnx package has what seem like excellent tools to tool to calculate "intersection density," but for specified polygons, like city boundaries, and does not seem to have functionality to create calls within Python by way or node ID.
I also know that I could probably do this in ArcGIS or QGIS, but because I have these OSM IDs handy in my database already, it just seems like a waste to load a whole dang shapefile of intersections for a large metropolitan region and do some kind of buffering process to get the information I need. Plus if I had a script handy to extract some pieces of information by way or node ID, I could more easily broaden the kinds of data I extract to get other tidbits of great information recorded for OSM elements.
Thank you spatial data community!
Traffic signals should always be tagged "highway" = "traffic_signals", but individual nodes may also be tagged with a key of "traffic_signals". Thus the first step to get all traffic signals can be done like this:
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]))
That then holds all the OSM ID values of nodes describing traffic signals.
One straightforward way to obtain junction densities is to convert the sf representation of highways to a dodgr network, which is a simple data.frame with each row being a network edge. The poly2line step converts strict sf polygons such as roundabouts into linestring objects, while the dodgr_contract_graph() call reduces the network to junction vertices only.
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) will give you the frequencies of different kinds of highways. You can then inspect junction frequencies for particular types as follows
net_footway <- net[net$highway == "footway", ]
table(table(c(net_footway$from_id, net_footway$to_id)))
Values of 1 indicate one-way terminal nodes; values of 2 indicate two-way terminal nodes; values of 4 indicate cross-junctions between two edges; and so on. That table goes up to 14, because footways can be quite complex, and there is obviously a junction of seven footways somewhere in Minneapolis. Those IDs are OSM ids, so you can easily check which of them are in the signal_node_ids values to determine which have traffic signals.
Remaining issues to address:
"highway" types are easy, but intersections between different types would require more complicated modifications to this code. Although straightforward, you'd need to consistently sub-set the dodgr data.frame in which edges are directed: $from_id -> $to_id."traffic_signals". Note that there is no "right" way to do this, because, for example, junctions may have separate signals for pedestrians and automobiles, and decisions of whether to consider these to be "the same" signals will always entail some degree of subjectivity.