Calculating shortest path from and to any spatial point

[luukvdmeer]

Is your feature request related to a problem? Please describe.
In tidygraph and igraph, you have to provide a node index as from and to in shortest path calculation. This requires you to know your node indices. In spatial networks, it is common to calculate the shortest path from and to a spatial point with X and Y coordinates, no matter if this point is a node in the network or not.

Describe the solution you'd like
Enable spatial shortest path calculation as described above. The function will find either the nearest node (with st_snap_point_to_node) or the nearest point on the nearest edge (with st_snap_point_to_edge followed by st_split_edge), and then calculate the shortest path as normal.

[luukvdmeer]

A first implementation of this is now available.

The function st_shortest_paths wraps shortest_paths from igraph:

library(sfnetworks)

net = as_sfnetwork(roxel, directed = FALSE)

st_shortest_paths(net, 1, 9)$vpath
#> [[1]]
#> + 17/701 vertices, from 1e53b31:
#>  [1]   1 644 579 260 164 208 426 424 533 657 512 506 111 524 292 167   9

igraph::shortest_paths(net, 1, 9)$vpath
#> [[1]]
#> + 17/701 vertices, from 1e53b31:
#>  [1]   1 644 579 260 164 208 426 424 533 657 512 506 111 524 292 167   9

What makes st_shortest_paths different, is that instead of providing node indices as from and to, you can provide sf or sfc objects.

p1 = net %>%
  activate(nodes) %>%
  st_as_sf() %>%
  slice(1)

p1
#>  Simple feature collection with 1 feature and 0 fields
#>  geometry type:  POINT
#>  dimension:      XY
#>  bbox:           xmin: 7.533722 ymin: 51.95556 xmax: 7.533722 ymax: 51.95556
#>  CRS:            EPSG:4326
#>  # A tibble: 1 x 1
#>                      x
#>            <POINT [°]>
#>  1 (7.533722 51.95556)

p2 = net %>%
  activate(nodes) %>%
  st_as_sf() %>%
  slice(9)

p2
#>  Simple feature collection with 1 feature and 0 fields
#>  geometry type:  POINT
#>  dimension:      XY
#>  bbox:           xmin: 7.537672 ymin: 51.9475 xmax: 7.537672 ymax: 51.9475
#>  CRS:            EPSG:4326
#>  # A tibble: 1 x 1
#>                     x
#>           <POINT [°]>
#>  1 (7.537673 51.9475)

st_shortest_paths(net, p1, p2)$vpath
#>  [[1]]
#>  + 17/701 vertices, from 1e53b31:
#>  [1]   1 644 579 260 164 208 426 424 533 657 512 506 111 524 292 167   9

In the case above, these points lie on the network. However, they don't have to. Points that don't lie on the network will be replaced by their nearest node.

lengths(sf::st_intersects(p1, net)) > 0
#> [1] TRUE

lengths(sf::st_intersects(p2, net)) > 0
#> [1] TRUE

p3 = sf::st_sfc(
  sf::st_geometry(p1)[[1]] + sf::st_point(c(0.01, 0.01)), 
  crs = sf::st_crs(p1)
)

p4 = sf::st_sfc(
  sf::st_geometry(p2)[[1]] + sf::st_point(c(0.01, 0.01)), 
  crs = sf::st_crs(p2)
)

p3
#> Geometry set for 1 feature 
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 7.543722 ymin: 51.96556 xmax: 7.543722 ymax: 51.96556
#> CRS:            EPSG:4326
#> POINT (7.543722 51.96556)

p4
#> Geometry set for 1 feature 
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 7.547672 ymin: 51.9575 xmax: 7.547672 ymax: 51.9575
#> CRS:            EPSG:4326
#> POINT (7.547673 51.9575)

lengths(sf::st_intersects(p3, net)) > 0
#> [1] FALSE

lengths(sf::st_intersects(p4, net)) > 0
#> [1] FALSE

st_shortest_paths(net, p3, p4)$vpath
#> [[1]]
#> + 6/701 vertices, from 1e53b31:
#> [1] 180 178 179 601 192 602

All paths above are calculated without edge weights. Just as in igraph::shortest_paths you can provide those weights through the weights argument. However, in st_shortest_path you also have the option to just provide a column name to the weights argument, which should match an existing column in the networks edges table, which will then be used as weights.

net  = net %>%
  activate(edges) %>%
  dplyr::mutate(length = sf::st_length(.))

net
#> # An sfnetwork with 701 nodes and 851 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # An undirected multigraph with 14 components and spatially explicit edges
#> #
#> # Edge Data:     851 x 6 (active)
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 7.522594 ymin: 51.94151 xmax: 7.546705 ymax:
#> #   51.9612
#>    from    to name      type                            geometry   length
#>   <int> <int> <fct>     <fct>                   <LINESTRING [°]>      [m]
#> 1     1     2 Havixbec… resid… (7.533722 51.95556, 7.533461 51.…  28.859…
#> 2     3     4 Pienersa… secon… (7.532442 51.95422, 7.53236 51.9… 107.704…
#> 3     5     6 Schulte-… resid… (7.532709 51.95209, 7.532823 51.…  54.362…
#> 4     7     8 NA        path   (7.540063 51.94468, 7.539696 51.… 155.230…
#> 5     9    10 Welsingh… resid… (7.537673 51.9475, 7.537614 51.9… 208.663…
#> 6    11    12 NA        footw… (7.543791 51.94733, 7.54369 51.9…  63.023…
#> # … with 845 more rows
#> #
#> # Node Data:     701 x 1
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 7.522622 ymin: 51.94151 xmax: 7.546705 ymax:
#> #   51.9612
#>                     x
#>           <POINT [°]>
#> 1 (7.533722 51.95556)
#> 2 (7.533461 51.95576)
#> 3 (7.532442 51.95422)
#> # … with 698 more rows

st_shortest_paths(net, p1, p2, weights = "length")$vpath
#> [[1]]
#> + 17/701 vertices, from 1e53b31:
#>  [1]   1 645 481 263 164 208 426 424 533 657 512 506 111 524 292 167   9

Alternatively, you can name your weights column weight. That is the same as having a weight attribute in an igraph object. The weight values will then be automatically detected and do not have to be provided explicitly.

net  = net %>%
  activate(edges) %>%
  dplyr::rename(weight = length)

st_shortest_paths(net, p1, p2)$vpath
#> [[1]]
#> + 17/701 vertices, from 1e53b31:
#>  [1]   1 645 481 263 164 208 426 424 533 657 512 506 111 524 292 167   9

I now showed 1-to-1 shortest paths. You can also calculate 1-to-many shortest paths, by giving a multiple feature sf or sfc object as to argument. Having multiple features as from is not possible in this case (since igraph does not support it).

Having multiple features as from is possible in st_network_distances, which is a wrapper around igraph::distances. This returns a distance matrix. I chose the name st_network_distances over st_distances to not confuse with the sf function st_distance.

ps1 = rbind(p1, sf::st_as_sf(p3))
ps2 = rbind(p2, sf::st_as_sf(p4))

ps1
#> Simple feature collection with 2 features and 0 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 7.533722 ymin: 51.95556 xmax: 7.543722 ymax: 51.96556
#> CRS:            EPSG:4326
#> # A tibble: 2 x 1
#>                     x
#>           <POINT [°]>
#> 1 (7.533722 51.95556)
#> 2 (7.543722 51.96556)

ps2
#> Simple feature collection with 2 features and 0 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 7.537672 ymin: 51.9475 xmax: 7.547672 ymax: 51.9575
#> CRS:            EPSG:4326
#> # A tibble: 2 x 1
#>                    x
#>          <POINT [°]>
#> 1 (7.537673 51.9475)
#> 2 (7.547673 51.9575)

# Note that our net object now has a weight column which will be automatically used.
st_network_distances(net, ps1, ps2)
#>          [,1]      [,2]
#> [1,] 1040.254 1074.5842
#> [2,] 1805.846  298.9539

Finally, there is also st_all_shortest_paths, which is a wrapper around igraph::all_shortest_paths. It works the same as st_shortest_paths, but in the case there are multiple shortest paths from A to B, it will return all of them, instead of only one.

Remark
If you provide points that are not network nodes, the nearest node to those points will be taken. I would like to also give the option to take the nearest point on the nearest edge. However, this seems to be not as simple as it sounds. I plan to open a hackathon issue about that.

[Robinlovelace]

This looks really great @luukvdmeer, excellent work. This links to a request made by @rafapereirabr back in 2017 to be able to add new nodes to the network for routing to the nearest point: ropensci/stplanr#237 (comment)

I think that will be a good hack and can be broken-up into modular parts:

• Identifying the point on the nearest edge that is closest to the point (@agila5 has looked at this I think, do you know of any implementation?)
• Adding new nodes to the network
• Routing to the new node
• Writing a function that combines the steps outlined above

[luukvdmeer]

I wrote a function that find the point on the nearest edge that is closest to the input point:

# x is a set of points.
# y is the edges element of the sfnetwork, containing lines.
st_snap_to_network = function(x, y) {
  # Function to find the nearest edge point to a single input point.
  f = function(z) {
    # Convert sfg to sfc.
    z = sf::st_sfc(z, crs = sf::st_crs(x))
    # Run st_nearest_points.
    # This returns for each combination (z, y) a linestring from z to the
    # .. nearest point on y.
    np = sf::st_nearest_points(z, y)
    # Then, the endpoint of the np-line to the nearest feature y to z
    # .. is the nearest edge point to z.
    lwgeom::st_endpoint(np[sf::st_nearest_feature(z, y),])
  }
  # Run f for all input points.
  geoms = do.call("c", lapply(sf::st_geometry(x), f))
  # Replace the geometries of x with the corresponding nearest edge points.
  if (inherits(x, "sf")) st_geometry(x) = geoms
  if (inherits(x, "sfc")) x = geoms
  x
}

However, there is a problem with this approach. The returned point does not always intersect with the network.. Sometimes it is located a few millimeters away from the edge. This has to do with floating point precision.

See https://gis.stackexchange.com/questions/288988/nearest-point-does-not-intersect-line-in-sf-r-package
And r-spatial/sf#790

This reply does not sound hopeful:

I think this problem cannot be solved. Setting precision essentially rounds coordinates to a (scaled) integer grid. If you want to represent points on lines connecting points on an integer grid exactly, you need to store them as rational numbers, i.e. as a ratio of two (possibly very large) integers. What R (or GEOS) does is approximating this ratio when storing it as an eight byte double. By that, we're back at R FAQ 7.31.

There is a recent blogpost addressing this problem as well, which uses a somehow "hacky" solution: https://ryanpeek.org/mapping-in-R-workshop/03_vig_snapping_points_to_line.html

A lot to think and test, ideal for the hackathon :)

[Robinlovelace]

Very nice, and easier than I thought to get the nearest point on a geometry. Just realised this after re-running great example in sf:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
r = sqrt(2)/10
pt1 = st_point(c(.1,.1))
pt2 = st_point(c(.9,.9))
pt3 = st_point(c(.9,.1))
b1 = st_buffer(pt1, r)
b2 = st_buffer(pt2, r)
b3 = st_buffer(pt3, r)
(ls0 = st_nearest_points(b1, b2)) # sfg
(ls = st_nearest_points(st_sfc(b1), st_sfc(b2, b3))) # sfc
plot(b1, xlim = c(-.2,1.2), ylim = c(-.2,1.2), col = NA, border = 'green')
plot(st_sfc(b2, b3), add = TRUE, col = NA, border = 'blue')
plot(ls, add = TRUE, col = 'red')

^{Created on 2020-06-08 by the reprex package (v0.3.0)}

This will be great for a hackathon and it looks like a hard part is already solved!

[rafapereirabr]

Thank you Robin for the tagging me here. Lucas, This is looking really good. Congrats on the work! I've registered for the webinar and I'm looking forward to seeing the progress with sfnetworks.

[luukvdmeer]

Pushed to master. For the nearest_point_on_edge discussion, see #54