Get shortest/fastest path between two points from SpatialLines(DataFrame)

A way to approach this problem could be by switching from sp to sf objects, and then use sfnetworks and tidygraph which internally use igraph.

First we recreate the lines and points as sf objects.

x <- c(1,5,6,8)
y1 <- c(1,3,4,7)
y2 <- c(5,5,5,2)

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
L = st_sf(
  ID = c("a", "b"),
  geometry = c(
    st_sfc(st_linestring(cbind(x,y1))),
    st_sfc(st_linestring(cbind(x,y2)))
  )
)
P = st_as_sf(data.frame(x=c(1,8),y=c(1,2)), coords = c("x","y"))
P_snap = st_as_sf(data.frame(x=8,y=1), coords = c("x","y"))

Then we can create an sfnetwork object from the LINESTRINGs L.

library(sfnetworks)
(G = L %>%
  as_sfnetwork(directed = F))
#> # A sfnetwork with 4 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # An unrooted forest with 2 trees with spatially explicit edges
#> #
#> # Node Data:     4 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 1 xmax: 8 ymax: 7
#>   geometry
#>    <POINT>
#> 1    (1 1)
#> 2    (8 7)
#> 3    (1 5)
#> 4    (8 2)
#> #
#> # Edge Data:     2 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 1 xmax: 8 ymax: 7
#>    from    to ID                geometry
#>   <int> <int> <chr>         <LINESTRING>
#> 1     1     2 a     (1 1, 5 3, 6 4, 8 7)
#> 2     3     4 b     (1 5, 5 5, 6 5, 8 2)
G %>% class()
#> [1] "sfnetwork" "tbl_graph" "igraph"
plot(G)

This will return a network consisting of the two linestrings as edges, and nodes created at the linestring end points. However, routing will be difficult since there is no node that connects both edges. This will probably not be a problem in the larger river network, but for this particular case we can fix ti with the following code:

G = L %>%
  # Combine LINESTRINGS into a MULTILINESTRING geometry
  st_combine() %>%
  # Create a node where the MULTILINESTRINGs cross
  st_node() %>%
  # Cast back to LINESTRINGs
  st_cast('LINESTRING') %>%
  # Reset IDs
  st_sf(ID = c("a","a","b","b")) %>%
  # Create sfnetwork
  as_sfnetwork(directed = F)

Now, with help of tidygraph we can convert the network to its spatial shortest paths. Convert will result in a subset of the original network, consisting only of the spatial shortest path in this case. sfnetworks accepts sf POINTs as to and from arguments, while internally snapping this points to the nearest node. It also uses the edge geographical length as weight internally

library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
# Convert to spatial shortest path between points in the network
G_sp1 = G %>%
  convert(to_spatial_shortest_paths, from = P[1,], to = P[2,])
# Convert to spatial shortest path to snapped point
G_sp2 = G %>%
  convert(to_spatial_shortest_paths, from = P[1,], to = P_snap)

We can plot the result to see how it worked

par(mar = c(0,0,0,0), mfrow = c(1,2))

plot(G, col = "grey", reset = FALSE)
plot(P, pch = 8, col = "green", add = TRUE)
G_sp1 %>%
  plot(col = "blue", add = TRUE)

plot(G, col = "grey", reset = FALSE)
plot(P, pch = 8, col = "green", add = TRUE)
plot(P_snap, pch = 8, col = "red", add = TRUE)
G_sp2 %>%
  plot(col = "purple", add = TRUE)

If you still need to obtain the shortest path as a SpatialLine, it can be extracted as follows:

(L_sp1 = G_sp1 %>%
  st_geometry(active = "edges") %>%
  as_Spatial())
#> An object of class "SpatialLines"
#> Slot "lines":
#> [[1]]
#> An object of class "Lines"
#> Slot "Lines":
#> [[1]]
#> An object of class "Line"
#> Slot "coords":
#>          [,1] [,2]
#> [1,] 1.000000  1.0
#> [2,] 5.000000  3.0
#> [3,] 6.000000  4.0
#> [4,] 6.333333  4.5
#> 
#> 
#> 
#> Slot "ID":
#> [1] "ID1"
#> 
#> 
#> [[2]]
#> An object of class "Lines"
#> Slot "Lines":
#> [[1]]
#> An object of class "Line"
#> Slot "coords":
#>          [,1] [,2]
#> [1,] 6.333333  4.5
#> [2,] 8.000000  2.0
#> 
#> 
#> 
#> Slot "ID":
#> [1] "ID2"
#> 
#> 
#> 
#> Slot "bbox":
#>   min max
#> x   1 8.0
#> y   1 4.5
#> 
#> Slot "proj4string":
#> CRS arguments: NA