1.1 Dual-Weighted Directed Graphs

Shortest-path distances on weighted graphs can be calculated using a number of other R packages, such as igraph or e1071. dodgr comes into its own through its ability to trace paths through dual-weighted directed graphs, illustrated in Fig. 2.

Dual-weighted directed graphs are common in many areas, a foremost example being routing through street networks. Routes through street networks depends on mode of transport: the route a pedestrian might take will generally differ markedly from the route the same person might take if behind the wheel of an automobile. Routing through street networks thus generally requires each edge to be specified with two weights or distances: one quantifying the physical distance, and a second weighted version reflecting the mode of transport (or some other preferential weighting).

dodgr calculates shortest paths using one set of weights (called “weights” or anything else starting with “w”), but returns the actual lengths of them using a second set of weights (called “distances”, or anything else starting with “d”). If no weights are specified, distances alone are used both for routing and final distance calculations. Consider that the weights and distances of Fig. 2 are the black and grey lines, respectively, with the latter all equal to one. In this case, the graph and associated shortest distances are,

##   from to w d
## 1    A  B 1 1
## 2    B  A 2 1
## 3    B  C 1 1
## 4    B  D 3 1
## 5    C  B 2 1
## 6    C  D 1 1
## 7    D  C 2 1
## 8    D  A 1 1

##   A B C D
## A 0 1 2 2
## B 1 0 1 1
## C 2 1 0 1
## D 1 2 1 0

Note that even though the shortest “distance” from A to D is actually A→B→D with a distance of only 2, that path has a weighted distance of 1 + 3 = 4. The shortest weighted path is A→B→C→D, with a distance both weighted and unweighted of 1 + 1 + 1 = 3. Thus d(A,D) = 3 and not 2.

3.1 Graph Components

The foregoing graph object returned from weight_streetnet() also includes a $component column enumerating all of the distinct inter-connected components of the graph.

table (graph$component)

## 
##    1    2    3 
## 4649 2066   98

Components are numbered in order of decreasing size, with $component = 1 always denoting the largest component. In this case, that component contains 3,934 edges, representing 65% of the graph. There are clearly only three distinct components, but this number may be much larger for larger graphs, and may be obtained from,

length (unique (graph$component))

## [1] 3

Component numbers can be determined for any types of graph with the dodgr_components() function. For example, the following lines reduce the previous graph to a minimal (non-spatial) structure of four columns, and then (re-)calculate a fifth column of $components:

cols <- c ("edge_id", "from_id", "to_id", "d")
graph_min <- graph [, which (names (graph) %in% cols)]
graph_min <- dodgr_components (graph_min)
head (graph_min)

##   edge_id    from_id      to_id          d component
## 1       1  339318500  339318502 130.000241         1
## 2       2  339318502  339318500 130.000241         1
## 3       3  339318502 2398958028   8.890622         1
## 4       4 2398958028  339318502   8.890622         1
## 5       5 2398958028 1427116077   9.307736         1
## 6       6 1427116077 2398958028   9.307736         1

The component column column can be used to select or filter any component in a graph. It is particularly useful to ensure routing calculations consider only connected vertices through simply removing all minor components:

graph_connected <- graph [graph$component == 1, ]

This is explored further below (under Distance Matrices).

3.2 Weighting Profiles

Dual-weights for street networks are generally obtained by multiplying the distance of each segment by a weighting factor reflecting the type of highway. As demonstrated above, this can be done easily within dodgr with the weight_streetnet() function, which applies the named weighting profiles included with the dodgr package to OpenStreetMap networks extracted with the osmdata package.

This function uses the internal data dodgr::weighting_profiles, which is a list of three items:

1. weighting_profiles;
2. surface_speeds; and
3. penalties

Most of these data are used to calculate routing times with the dodgr_times function, as detailed in an additional vignette. The only aspects relevant for distances are the profiles themselves, which assign preferential weights to each distinct type of highway.

wp <- weighting_profiles$weighting_profiles
names (wp)

## [1] "name"      "way"       "value"     "max_speed"

class (wp)

## [1] "data.frame"

unique (wp$name)

##  [1] "foot"       "horse"      "wheelchair" "bicycle"    "moped"     
##  [6] "motorcycle" "motorcar"   "goods"      "hgv"        "psv"

wp [wp$name == "foot", ]

##    name            way value max_speed
## 1  foot       motorway  0.00        NA
## 2  foot          trunk  0.40        NA
## 3  foot        primary  0.50         5
## 4  foot      secondary  0.60         5
## 5  foot       tertiary  0.70         5
## 6  foot   unclassified  0.80         5
## 7  foot    residential  0.90         5
## 8  foot        service  0.90         5
## 9  foot          track  0.95         5
## 10 foot       cycleway  0.95         5
## 11 foot           path  1.00         5
## 12 foot          steps  0.80         2
## 13 foot          ferry  0.20         5
## 14 foot  living_street  0.95         5
## 15 foot      bridleway  1.00         5
## 16 foot        footway  1.00         5
## 17 foot     pedestrian  1.00         5
## 18 foot  motorway_link  0.00        NA
## 19 foot     trunk_link  0.40        NA
## 20 foot   primary_link  0.50         5
## 21 foot secondary_link  0.60         5
## 22 foot  tertiary_link  0.70         5

Each profile is defined by a series of percentage weights quantifying highway-type preferences for a particular mode of travel. The distinct types of highways within the Hampi graph obtained above can be tabulated with:

table (graph$highway)

## 
## living_street          path       primary   residential     secondary 
##            20          3557           430           196           560 
##       service         steps         track  unclassified 
##           256           108           914           772

Hampi is unlike most other human settlements on the planet in being a Unesco World Heritage area in which automobiles are generally prohibited. Accordingly, numbers of "footway", "path", and "pedestrian" ways far exceed typical categories denoting automobile traffic ("primary", "residential", "tertiary")

It is also possible to use other types of (non-OpenStreetMap) street networks, an example of which is the os_roads_bristol data provided with the package. “OS” is the U.K. Ordnance Survey, and these data are provided as a Simple Features (sf) data.frame with a decidedly different structure to osmdata data.frame objects:

names (hampi) # many fields manually removed to reduce size of this object

##  [1] "osm_id"        "bicycle"       "covered"       "foot"         
##  [5] "highway"       "incline"       "motorcar"      "motorcycle"   
##  [9] "motor_vehicle" "oneway"        "surface"       "tracktype"    
## [13] "tunnel"        "width"         "geometry"

names (os_roads_bristol)

##  [1] "fictitious" "identifier" "class"      "roadNumber" "name1"     
##  [6] "name1_lang" "name2"      "name2_lang" "formOfWay"  "length"    
## [11] "primary"    "trunkRoad"  "loop"       "startNode"  "endNode"   
## [16] "structure"  "nameTOID"   "numberTOID" "function."  "geometry"

The latter may be converted to a dodgr network by first specifying a weighting profile, here based on the formOfWay column:

colnm <- "formOfWay"
table (os_roads_bristol [[colnm]])

## 
## Collapsed Dual Carriageway           Dual Carriageway 
##                         14                          6 
##         Single Carriageway                  Slip Road 
##                          1                          8

wts <- data.frame (name = "custom",
                   way = unique (os_roads_bristol [[colnm]]),
                   value = c (0.1, 0.2, 0.8, 1))
net <- weight_streetnet (os_roads_bristol, wt_profile = wts,
                         type_col = colnm, id_col = "identifier")

The resultant net object contains the street network of os_roads_bristol weighted by the specified profile, and in a format suitable for submission to any dodgr routine.

3.3 Random Sub-Graphs

The dodgr packages includes a function to select a random connected portion of graph including a specified number of vertices. This function is used in the compare_heaps() function described below, but is also useful for general statistical analyses of large graphs which may otherwise take too long to compute.

graph_sub <- dodgr_sample (graph, nverts = 100)
nrow (graph_sub)

## [1] 199

The random sample has around twice as many edges as vertices, in accordance with the statistics calculated above.

nrow (dodgr_vertices (graph_sub))

## [1] 100

4.1 Aligning Routing Points to Graphs

For spatial graphs—that is, those containing columns of latitudes and longitudes (or “x” and “y”)—routing points can be represented by a simple matrix of arbitrary latitudes and longitudes (or, again, “x” and “y”). dodgr_distances() will map these points to the closest network points, and return corresponding shortest-path distances. This may be illustrated by generating random points within the bounding box of the above map of Hampi. As demonstrated above, the coordinates of all vertices may be extracted with the dodgr_vertices() function, enabling random points to be generated with the following lines:

vt <- dodgr_vertices (graph)
n <- 100 # number of points to generate
xy <- data.frame (x = min (vt$x) + runif (n) * diff (range (vt$x)),
                  y = min (vt$y) + runif (n) * diff (range (vt$y)))

Submitting these to dodgr_distances() as points from which to route will generate a distance matrix from each of these 100 points to every other point in the graph:

d <- dodgr_dists (graph, from = xy)
dim (d); range (d, na.rm = TRUE)

## [1]  100 3337

## [1]     0.00 14926.04

If the to argument is also specified, the matrix returned will have rows matching from and columns matching to

d <- dodgr_dists (graph, from = xy, to = xy [1:10, ])
dim (d)

## [1] 100  10

Some of the resultant distances in the above cases are NA because the points were sampled from the entire bounding box, and the street network near the boundaries may be cut off from the rest. As demonstrated above, the weight_streetnet() function returns a component vector, and such disconnected edges will have graph$component > 1, because graph$component == 1 always denotes the largest connected component. This means that the graph can always be reduced to the single largest component with the following single line:

graph_connected <- graph [graph$component == 1, ]

A distance matrix obtained from running dodgr_distances on graph_connected should generally contain no NA values, although some points may still be effectively unreachable due to one-way connections (or streets). Thus, routing on the largest connected component of a directed graph ought to be expected to yield the minimal number of NA values, which may sometimes be more than zero. Note further that spatial routing points (expressed as from and/or to arguments) will in this case be mapped to the nearest vertices of graph_connected, rather than the potentially closer nearest points of the full graph. This may make the spatial mapping of routing points less accurate than results obtained by repeating extraction of the street network using an expanded bounding box. For automatic extraction of street networks with dodgr_distances(), the extent by which the bounding box exceeds the range of routing points (from and to arguments) is determined by an extra parameter expand, quantifying the relative extent to which the bounding box should exceed the spatial range of the routing points. This is illustrated in the following code which calculates distances between 100 random points:

bb <- osmdata::getbb ("york uk")
npts <- 100
xy <- apply (bb, 1, function (i) min (i) + runif (npts) * diff (i))

routed_points <- function (expand = 0, pts) {

    gr0 <- dodgr_streetnet (pts = pts, expand = expand) %>%
        weight_streetnet ()
    d0 <- dodgr_dists (gr0, from = pts)
    length (which (is.na (d0))) / length (d0)
}

vapply (c (0, 0.05, 0.1, 0.2), function (i) routed_points (i, pts = xy),
        numeric (1))

## [1] 0.04007477 0.02326452 0.02131992 0.00000000

With a street network that precisely encompasses the submitted routing points (expand = 0), 4% of pairwise distances are unable to be calculated; with a bounding box expanded to 5% larger than the submitted points, this is reduced to 2.3%, and with expansion to 20%, all points can be connected.

For non-spatial graphs, from and to must match precisely on to vertices named in the graph itself. In the graph considered above, these vertex names were contained in the columns, from_id and to_id. The minimum that a dodgr graph requires is,

head (graph [, names (graph) %in% c ("from_id", "to_id", "d")])

##      from_id      to_id          d
## 1  339318500  339318502 130.000241
## 2  339318502  339318500 130.000241
## 3  339318502 2398958028   8.890622
## 4 2398958028  339318502   8.890622
## 5 2398958028 1427116077   9.307736
## 6 1427116077 2398958028   9.307736

in which case the from values submitted to dodgr_dists() (and to, if given) must directly name the vertices in the from_id and to_id columns of the graph. This is illustrated in the following code:

graph_min <- graph [, names (graph) %in% c ("from_id", "to_id", "d")]
fr <- sample (graph_min$from_id, size = 10) # 10 random points
to <- sample (graph_min$to_id, size = 20)
d <- dodgr_dists (graph_min, from = fr, to = to)
dim (d)

## [1] 10 20

The result is a 10-by-20 matrix of distances between these named graph vertices.

4.2 Shortest Path Calculations: Priority Queues

dodgr uses an internal library Shane Saunders (2004) for the calculation of shortest paths using a variety of priority queues (see Miller 1960 for an overview). In the context of shortest paths, priority queues determine the order in which a graph is traversed (Tarjan 1983), and the choice of priority queue can have a considerable effect on computational efficiency for different kinds of graphs (Johnson 1977). In contrast to dodgr, most other R packages for shortest path calculations do not use priority queues, and so may often be less efficient. Shortest path distances can be calculated in dodgr with priority queues that use the following heaps:

1. Binary heaps;
2. Fibonacci heaps (Fredman and Tarjan 1987);
3. Trinomial and extended trinomial heaps (Takaoka 2000); and
4. 2-3 heaps (Takaoka 1999).

Differences in how these heaps operate are often largely extraneous to direct application of routing algorithms, even though heap choice may strongly affect performance. To avoid users needing to know anything about algorithmic details, dodgr provides a function compare_heaps() to which a particular graph may be submitted in order to determine the optimal kind of heap.

The comparisons are actually made on a randomly selected sub-component of the graph containing a defined number of vertices (with a default of 1,000, or the entire graph if it contains fewer than 1,000 vertices).

compare_heaps (graph, nverts = 100)

## Loading required namespace: igraph

## # A tibble: 11 × 6
##    expression                 min   median `itr/sec` mem_alloc `gc/sec`
##    <bch:expr>            <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
##  1 BHeap                   1.47ms   1.52ms      655.   35.83KB     10.5
##  2 FHeap                   1.48ms   1.54ms      637.   35.83KB     10.6
##  3 TriHeap                 1.48ms   1.53ms      651.   35.83KB     10.7
##  4 TriHeapExt              1.21ms   1.26ms      792.   38.87KB     10.5
##  5 Heap23                  1.48ms   1.55ms      640.   35.83KB     10.6
##  6 BHeap_contracted        1.26ms   1.31ms      752.    5.44KB     12.7
##  7 FHeap_contracted        1.25ms   1.33ms      749.    5.44KB     10.6
##  8 TriHeap_contracted      1.27ms   1.34ms      735.    5.44KB     10.6
##  9 TriHeapExt_contracted 961.45µs   1.02ms      933.    5.44KB     15.1
## 10 Heap23_contracted       1.26ms   1.32ms      752.    5.44KB     10.5
## 11 igraph                621.22µs 662.57µs     1465.  453.52KB     16.0

The key column of that data.frame is relative, which quantifies the relative performance of each test in relation to the best which is given a score of 1. dodgr using the default heap = "BHeap", which is a binary heap priority queue, performs faster than igraph (Csardi and Nepusz 2006) for these graphs. Different kind of graphs will perform differently with different priority queue structures, and this function enables users to empirically discern the optimal heap for their kind of graph.

Note, however, that this is not an entirely fair comparison, because dodgr calculates dual-weighted distances, whereas igraph—and indeed all other R packages—only directly calculate distances based on a single set of weights. Implementing dual-weighted routing in these cases requires explicitly re-tracing all paths and summing the second set of weights along each path. A time comparison in that case would be very strongly in favour of dodgr. Moreover, dodgr can convert graphs to contracted form through removing redundant vertices, as detailed in the following section. Doing so greatly improves performance with respect to igraph.

For those wishing to do explicit comparisons themselves, the following code generates the igraph equivalent of dodgr_distances(), although of course for single-weighted graphs only:

v <- dodgr_vertices (graph)
pts <- sample (v$id, 1000)
igr <- dodgr_to_igraph (graph)
d <- igraph::distances (igr, v = pts, to = pts, mode = "out")

5.1 Routing on Contracted Graphs

Routing is often desired between defined points, and these points may inadvertently be removed in graph contraction. The dodgr_contract_graph() function accepts an additional argument specifying vertices to keep within the contracted graph. This list of vertices must directly match the vertex ID values in the graph.

The following code illustrates how to retain specific vertices within contracted graphs:

grc <- dodgr_contract_graph (graph)
nrow (grc)

## [1] 748

verts <- sample (dodgr_vertices (graph)$id, size = 100)
head (verts) # a character vector

## [1] "2588146032" "2398957864" "2627422261" "7799710735" "6161853534"
## [6] "5351820914"

grc <- dodgr_contract_graph (graph, verts)
nrow (grc)

## [1] 941

Retaining the nominated vertices yields a graph with considerably more edges than the fully contracted graph excluding these vertices. The dodgr_distances() function can be applied to the latter graph to obtain accurate distances precisely routed between these points, yet using the speed advantages of graph contraction.