Network pre-processing instructions

[Robinlovelace]

After conversations in the Geocomputation with R issue tracker, and thanks to your input @vlarmet, I'm keen to give the package a try. However, I have no idea how to get started, beginning with just an OSM dataset, e.g. downloaded and imported into R with osmextract.

Here's a minimal reproducible example:

library(osmextract)
osm_lines = oe_get("Isle of Wight", stringsAsFactors = FALSE, quiet = TRUE)
osm_points = oe_get("Isle of Wight", layer = "points", stringsAsFactors = FALSE, quiet = TRUE)
nrow(osm_lines)
#> [1] 48281
nrow(osm_points)
#> [1] 61329
par(mar = rep(0, 4))
plot(st_geometry(osm_lines), xlim = c(-1.59, -1.1), ylim = c(50.5, 50.8))
plot(st_geometry(osm_points), xlim = c(-1.59, -1.1), ylim = c(50.5, 50.8))

I'm keen to contribute to this package if that would be helpful, could well be useful for my work...

[Robinlovelace]

Demo of this below. Are there any helper functions, e.g. to convert to/from sf?

library(dodgr)
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
library(tidyverse)
library(cppRouting)

hampi
#> Simple feature collection with 236 features and 14 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 76.37261 ymin: 15.30104 xmax: 76.49232 ymax: 15.36033
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>            osm_id bicycle covered foot    highway incline motorcar motorcycle
#> 28565950 28565950    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30643853 30643853     yes    <NA>  yes       path    <NA>      yes       <NA>
#> 30643907 30643907    <NA>    <NA>  yes       path    <NA>     <NA>       <NA>
#> 30659140 30659140    <NA>    <NA> <NA>    service    <NA>     <NA>       <NA>
#> 30659141 30659141    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30659142 30659142    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30659143 30659143    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30659150 30659150    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30663027 30663027    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#> 30663412 30663412    <NA>    <NA> <NA> pedestrian    <NA>     <NA>       <NA>
#>          motor_vehicle oneway     surface tracktype tunnel width
#> 28565950       private   <NA>        <NA>      <NA>   <NA>  <NA>
#> 30643853          <NA>   <NA>        sand      <NA>   <NA>  <NA>
#> 30643907          <NA>   <NA>        <NA>      <NA>   <NA>  <NA>
#> 30659140          <NA>   <NA>     asphalt      <NA>   <NA>  <NA>
#> 30659141          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#> 30659142          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#> 30659143          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#> 30659150          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#> 30663027          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#> 30663412          <NA>   <NA> cobblestone      <NA>   <NA>  <NA>
#>                                geometry
#> 28565950 LINESTRING (76.47491 15.341...
#> 30643853 LINESTRING (76.47398 15.312...
#> 30643907 LINESTRING (76.46027 15.330...
#> 30659140 LINESTRING (76.46524 15.334...
#> 30659141 LINESTRING (76.46116 15.334...
#> 30659142 LINESTRING (76.46091 15.336...
#> 30659143 LINESTRING (76.46062 15.335...
#> 30659150 LINESTRING (76.4609 15.3356...
#> 30663027 LINESTRING (76.46194 15.336...
#> 30663412 LINESTRING (76.46019 15.335...
class(hampi)
#> [1] "sf"         "data.frame"
plot(hampi)
#> Warning: plotting the first 9 out of 14 attributes; use max.plot = 14 to plot
#> all

graph = weight_streetnet (hampi, wt_profile = "foot")
head(graph)
#>   geom_num edge_id    from_id from_lon from_lat      to_id   to_lon   to_lat
#> 1        1       1  339318500 76.47491 15.34167  339318502 76.47612 15.34173
#> 2        1       2  339318502 76.47612 15.34173  339318500 76.47491 15.34167
#> 3        1       3  339318502 76.47612 15.34173 2398958028 76.47621 15.34174
#> 4        1       4 2398958028 76.47621 15.34174  339318502 76.47612 15.34173
#> 5        1       5 2398958028 76.47621 15.34174 1427116077 76.47628 15.34179
#> 6        1       6 1427116077 76.47628 15.34179 2398958028 76.47621 15.34174
#>            d d_weighted highway   way_id component      time time_weighted
#> 1 130.000241 130.000241    path 28565950         1 93.600174     93.600174
#> 2 130.000241 130.000241    path 28565950         1 93.600174     93.600174
#> 3   8.890622   8.890622    path 28565950         1  6.401248      6.401248
#> 4   8.890622   8.890622    path 28565950         1  6.401248      6.401248
#> 5   9.307736   9.307736    path 28565950         1  6.701570      6.701570
#> 6   9.307736   9.307736    path 28565950         1  6.701570      6.701570
v = dodgr_vertices(graph)
head(v)
#>            id        x        y component n
#> 1   339318500 76.47491 15.34167         1 0
#> 2   339318502 76.47612 15.34173         1 1
#> 4  2398958028 76.47621 15.34174         1 2
#> 6  1427116077 76.47628 15.34179         1 3
#> 8  7799710916 76.47634 15.34184         1 4
#> 10  339318503 76.47641 15.34190         1 5

graph_df = graph
class(graph_df) = "data.frame"
class(graph_df)
#> [1] "data.frame"

roads = graph_df |>
  filter(component == 1) |>
  transmute(from = from_id, to = to_id, weight = d_weighted)
coord = v |>
  transmute(ID = id, X = x, Y = y)

graph  =  makegraph(roads, directed = T, coords = coord)
graph$nbnode
#> [1] 2270

# from_id = roads$ID[1]
# to_id = roads$ID[2]
# path_pair1 = get_path_pair(graph, from = from_id, to = to_id)
origin  <-  sample(graph$dict$ref,  2)
destination  <-  sample(graph$dict$ref,  2)

ds = get_distance_pair(graph, origin, destination, algorithm = "NBA", constant = 110/0.06)
paths = get_path_pair(graph, origin, destination)
str(paths[[1]])
#>  chr [1:80] "2398957829" "338904628" "2398957830" "1054866750" "2398957833" ...
names(graph)
#> [1] "data"   "coords" "nbnode" "dict"   "attrib"
summary(paths[[1]] %in% graph$coords$ID)
#>    Mode    TRUE 
#> logical      80
path1_matrix = graph$coords[match(paths[[1]], graph$coords$ID), ]
path1_sf = sfheaders::sf_linestring(path1_matrix[, 2:3])
plot(path1_sf)

^{Created on 2023-07-01 with reprex v2.0.2}

[vlarmet]

Well, sfnetworks can do this right ?

[Robinlovelace]

Yes, but can it do it efficiently? Not sure. I think the dodgr approach of going from sf objects to the data frame format that cppRouting needs works well. I'm interested in what you'd recommend. I see you also recommend https://osm2po.de/ in your README. What would your recommended OSM -> routable network process be?

Happy to provide documentation on that and possibly even reproducible code (even if not in R).

[vlarmet]

Ok I see. To be honnest, I think it's an essential feature and I want to develop it from the beginning of cppRouting. cppRouting is a generic package so it should not be OSM specific. It's up to the user to update edge weights according to weighting profiles.

The main steps are :

• find starting/ending vertices of each LINESTRING
• create a set of (unique) nodes, based on coordinates
• match each edge starting/ending nodes (I think using distance through RANN package, exact matching seems to be risky for inaccurate shapefiles)
• function arguments must include : directed/undirected attribute for each edge, columns to keep (weights, length, anything), optionnal edge ID (if not, an integer ID is created)
• return network topology + edge ID + columns to keep

At first, input should be an sf object. In my dreams, cppRouting could extract network from a large .pbf file just like osm2po but it's a huge work.

[Robinlovelace]

Thanks Vincent. Big task as you say, but doable and important. Sounds like osm2po does a lot of the work, are you using that to produce the input graphs? Lots to discuss so may be worth chatting at some point in a quick video call? Will aim to send an email..

[ischlo]

Hi both, I might be a bit late on this one, but the best workflow using OSM networks and cppRouting I could come up with is using osmnx to download the data and then use the variables from, to and length from the edges layer together with the nodes layer to construct the graph. it's seamless and if you use reticulate for the python script, then you keep it all centralised in your R environment. osmnx is very good for the preprocessing of the OSM features into a 'graph like' structure required for cppRouting. Best, Ivann

[Robinlovelace]

Hey Ivann, thanks for sharing your experience. Easy to package-up/document? Wondering if a {cppRoutingsetup} type thing could be useful, certainly maybe for my use case although we're not using it for routing at present.

[ischlo]

Hi Robin yes, although I only use it locally for now. Indeed, it could be good to centralise a set of functions for the various use cases of cppR. Best, Ivann