R: sfnetworks: How to find routes between multiple A and B locations in a same dataset

Question

This is my dataset

df<-tribble(
  ~"shop.x",~"shop.y", ~"cust.x", ~"cust.y",
  78.100378,    9.944226,   78.096318,  9.954789,
  78.101155,    9.932190,   78.089824,  9.929975,
  78.141887,    9.928319,   78.110863,  9.952235,
  78.100381,    9.944226,   78.104066,  9.97013,
  78.097206,    9.948872,   78.11631,   9.947862
)

The df dataset has locations of shops and the customers.

I want to create shortest path route for each row (for every shop to customer location) using OSM map in R. Is it possible using sfnetworks?

The local road networks data is here

Answer 1

First, since this dataframe contains mainly coordinates in X and Y for shops and customers, I convert them into sf objects. Two objects for each set.

# Load packages
library(tidyverse)
library(sf)
library(sfnetworks)
library(tidygraph)

# Load the data
df <- tribble(
  ~"shop.x",~"shop.y", ~"cust.x", ~"cust.y",
  78.100378,    9.944226,   78.096318,  9.954789,
  78.101155,    9.932190,   78.089824,  9.929975,
  78.141887,    9.928319,   78.110863,  9.952235,
  78.100381,    9.944226,   78.104066,  9.97013,
  78.097206,    9.948872,   78.11631,   9.947862
)

# Convert into sf objects
shop = st_as_sf(df[,1:2], crs = 4326, coords = c("shop.x", "shop.y"))
cust = st_as_sf(df[,3:4], crs = 4326, coords = c("cust.x", "cust.y"))

The roads are LINESTRINGS which can be converted into a sfnetwork object. The network needs to be cleaned. More info on how to do network pre-processing here.

## I downloaded the data attached and saved it locally 
roads_dir = "data/roads.shp"
roads = st_read(roads_dir, crs = 4326)
#> Reading layer `Roads' from data source 
#> Simple feature collection with 6915 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 78.0217 ymin: 9.832762 xmax: 78.22482 ymax: 10.04193
#> Geodetic CRS:  WGS 84

# Convert into an sfnetwork
roads_sfn = as_sfnetwork(roads, directed = FALSE)

# Do some network cleaning
roads_clean = roads_sfn %>% 
  convert(to_spatial_subdivision, .clean = TRUE) %>% 
  convert(to_spatial_smooth, .clean = TRUE) 

To calculate the paths we can use an mapply function to the st_network_paths() function passing the from argument with the locations of the shops, and the to argument with the locations of the customers. Hence, the paths are calculated in order for each combination.

# Calculate the specified paths.
paths = mapply(
  st_network_paths,
  from = shop,
  to = cust,
  MoreArgs = list(x = roads_clean)
)["node_paths", ] %>%
  unlist(recursive = FALSE)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> Warning: Although argument from has length > 1, only the first element is used

The output is a list with five elements containing the node indices for each route, which correspond to the nodes in the original network.

paths[[1]]
#>  [1]  1842  1843 12191 12190   232  7977  7843 12174  8244  8205  7914  8139
#> [13]  8209  7958  7957  8229  8230  8328

If we want to get the path as a network itself, we can slice the original network as:

roads_clean %>%
  activate("nodes") %>%
  slice(paths[[1]])
#> # A sfnetwork with 18 nodes and 17 edges
#> #
#> # CRS:  EPSG:4326 
#> #
#> # An unrooted tree with spatially explicit edges
#> #
#> # Node Data:     18 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 78.09568 ymin: 9.943727 xmax: 78.10051 ymax: 9.955158
#>              geometry
#>           <POINT [°]>
#> 1 (78.10001 9.946228)
#> 2 (78.10051 9.943727)
#> 3 (78.10017 9.945039)
#> 4 (78.09918 9.951894)
#> 5 (78.09895 9.953372)
#> 6 (78.09647 9.955099)
#> # ... with 12 more rows
#> #
#> # Edge Data:     17 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 78.09568 ymin: 9.943727 xmax: 78.10051 ymax: 9.955158
#>    from    to                                                  geometry
#>   <int> <int>                                          <LINESTRING [°]>
#> 1     2     3 (78.10051 9.943727, 78.10029 9.944423, 78.10017 9.945039)
#> 2     6     7                    (78.09647 9.955099, 78.09877 9.954435)
#> 3    12    13                    (78.09569 9.955158, 78.09568 9.954763)
#> # ... with 14 more rows