How to convert a list of sf spatial points into a routable graph
I have an sf dataframe object with a series of points representing the shape of a bus route. I would like to turn this object into a routable graph so I can estimate the time it takes to traverse from point c to t.
Here is what I've tried using the dodgr package but I am not sure what I'm doing wrong here:
library(dodgr)
graph <- weight_streetnet(mydata, wt_profile = "motorcar", type_col="highway" , id_col = "id")
Error in check_highway_osmid(x, wt_profile) : Please specify type_col to be used for weighting streetnet
Reproducible data
The data looks like in the image below
mydata <- structure(list(shape_id = c(52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L,
52421L, 52421L, 52421L, 52421L, 52421L, 52421L, 52421L), length = structure(c(0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197, 0.191422504106197,
0.191422504106197, 0.191422504106197, 0.191422504106197), units = structure(list(
numerator = "km", denominator = character(0)), class = "symbolic_units"), class = "units"),
geometry = structure(list(structure(c(-46.5623281998182,
-23.5213458001468), class = c("XY", "POINT", "sfg")), structure(c(-46.562221,
-23.52129), class = c("XY", "POINT", "sfg")), structure(c(-46.562121,
-23.521235), class = c("XY", "POINT", "sfg")), structure(c(-46.5620233332577,
-23.5211840000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561925666591,
-23.5211330000609), class = c("XY", "POINT", "sfg")), structure(c(-46.561828,
-23.521082), class = c("XY", "POINT", "sfg")), structure(c(-46.5618098335317,
-23.5212126666783), class = c("XY", "POINT", "sfg")), structure(c(-46.5617916670273,
-23.5213433333544), class = c("XY", "POINT", "sfg")), structure(c(-46.5617735004869,
-23.5214740000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5617553339104,
-23.5216046667004), class = c("XY", "POINT", "sfg")), structure(c(-46.5617371672978,
-23.5217353333702), class = c("XY", "POINT", "sfg")), structure(c(-46.5617190006492,
-23.5218660000379), class = c("XY", "POINT", "sfg")), structure(c(-46.5617008339645,
-23.5219966667036), class = c("XY", "POINT", "sfg")), structure(c(-46.5616826672438,
-23.5221273333671), class = c("XY", "POINT", "sfg")), structure(c(-46.5616645004869,
-23.5222580000284), class = c("XY", "POINT", "sfg")), structure(c(-46.5616463336941,
-23.5223886666877), class = c("XY", "POINT", "sfg")), structure(c(-46.5616281668651,
-23.5225193333449), class = c("XY", "POINT", "sfg")), structure(c(-46.56161,
-23.52265), class = c("XY", "POINT", "sfg")), structure(c(-46.5617355000207,
-23.5226427501509), class = c("XY", "POINT", "sfg")), structure(c(-46.5618610000276,
-23.5226355002012), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = -46.5623281998182,
ymin = -23.52265, xmax = -46.56161, ymax = -23.521082), class = "bbox"), crs = structure(list(
epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L),
id = c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j",
"k", "l", "m", "n", "o", "p", "q", "r", "s", "t"), speed_kmh = c(11,
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
11, 11, 11, 11)), sf_column = "geometry", agr = structure(c(shape_id = NA_integer_,
length = NA_integer_, id = NA_integer_, speed_kmh = NA_integer_
), class = "factor", .Label = c("constant", "aggregate", "identity"
)), row.names = c("1.13", "1.14", "1.15", "1.16", "1.17", "1.18",
"1.19", "1.20", "1.21", "1.22", "1.23", "1.24", "1.25", "1.26",
"1.27", "1.28", "1.29", "1.30", "1.31", "1.32"), class = c("sf",
"data.table", "data.frame"))
The weight_streetnet function is really only designed to handle actual street networks, generally as produced by the osmdata::osmdata_sf/sp/sc() functions. It can nevertheless be tweaked to handle cases like this. The main thing needed is to convert the points into something that knows about edges in between them, like an sf::LINESTRING object:
x <- sf::st_combine (mydata) %>%
sf::st_cast ("LINESTRING") %>%
sf::st_sf ()
That gives a single-row object which can then be converted to dodgr format, and the id values matched back on to the edges
net <- weight_streetnet (x, type_col = "shape_id", id_col = "id", wt_profile = 1)
net$from_id <- mydata$id [as.integer (net$from_id)]
net$to_id <- mydata$id [as.integer (net$to_id)]
At that point, dodgr will have calculated and inserted the distances directly from the geographic coordinates. Your distances can then also be inserted and used for routing by replacing the d_weighted values:
net$d_weighted <- as.numeric (mydata$length [1])
dodgr_dists (net, from = "c", to = "t") # 236.0481
If you really want your distances to represent absolute distances used to calculate the final result, then simply replace the $d values
net$d <- net$d_weighted
dodgr_dists (net, from = "c", to = "t") # 3.254183
Note that for "simple" problems like this, igraph will generally be faster, because it calculates routes using a single set of weights. The only real advantage of dodgr in this context is the ability to use "dual weights" - the $d_weighted and $d values - such that the route is calculated according to $d_weighted, and the final distances according to $d.