How to show multiple roads using leaflet in R

Question

I am working on program connected to Traveling Salesman Problem. I use R and I load data where is lon,lng and calculate distances using osrm package. Like this:

distance_matrix <- matrix(NA, nrow = nrow(data), ncol = nrow(data))
for (i in 1:nrow(data)) {
  for (j in i:nrow(data)) {
    if (i != j) {
      route <- osrmRoute(c(data[i,]$lng,data[i,]$lat),c(data[j,]$lng,data[j,]$lat))
      distance_matrix[i,j] <- route$distance
      distance_matrix[j,i] <- route$distance
    }
  }
}

My question is- Is there any way after I calculate the TSP how to show final route connecting already calculated routes(and how to save them)?

I know how to show one route like this:

leaflet(route) |> 
  addTiles()  |>
  addPolylines(label = ~ sprintf("%d min - %d km", round(duration), round(distance)),
               labelOptions = labelOptions(noHide = FALSE))

but this show only one route and I would love to show after using the TSP algorithm the final smallest route(using already calculated one).

Thank u

Answer 1

If your aim is just to use a distance matrix and test out different TSP approaches, I'd guess it's fine to grab that matrix from osrm::osrmTable(), solve your task (assuming your solvers are returning a sequence of point names or indices, i.e. row numbers of distance matrix) and get a route from osrm::osrmRoute() for a location sequence. With some example data from osrm and TSP package it might look something like this:

library(osrm)
#> Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
#> Routing: OSRM - http://project-osrm.org/
library(leaflet)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

apotheke.df <- read.csv(system.file("csv/apotheke.csv", package = "osrm"))[1:10,]
apotheke.df
#>            id      lon      lat
#> 1   440338666 13.43853 52.47728
#> 2   538057637 13.57874 52.45823
#> 3   977657079 13.49039 52.61719
#> 4  3770254015 13.51974 52.49737
#> 5   364363337 13.45582 52.50112
#> 6  3666540067 13.33844 52.49479
#> 7  4234162689 13.35448 52.54872
#> 8  1492336686 13.43178 52.53216
#> 9  1298511709 13.38544 52.53797
#> 10  268726812 13.44518 52.47895

# get distnace matrices for duration and distance
distA <- osrmTable(loc = apotheke.df[c("lon", "lat")], measure = c('duration', 'distance'))

distA$durations
#>       1    2    3    4    5    6    7    8    9   10
#> 1   0.0 22.6 34.7 20.7 11.4 15.0 23.6 16.0 19.7  3.9
#> 2  23.2  0.0 42.4 15.8 20.0 31.8 38.6 26.6 33.2 20.3
#> 3  34.2 40.9  0.0 30.9 27.8 34.2 24.8 20.9 24.6 35.2
#> 4  20.4 15.2 30.1  0.0 12.9 28.3 28.1 16.1 22.8 18.4
#> 5  10.2 20.2 27.3 12.2  0.0 18.1 20.0  9.3 15.5  8.2
#> 6  14.2 31.5 34.1 28.2 18.3  0.0 14.6 16.7 14.6 15.4
#> 7  20.9 38.5 24.2 26.8 18.8 14.4  0.0 13.9  6.3 21.8
#> 8  15.1 27.5 20.3 15.7  9.3 16.4 14.7  0.0  9.0 16.1
#> 9  19.1 34.2 24.3 22.5 15.2 14.3  8.3  9.2  0.0 20.0
#> 10  4.3 20.3 34.8 17.6  8.3 17.0 25.1 16.1 21.0  0.0

distA$distances
#>        1     2     3     4     5     6     7     8     9    10
#> 1      0 13409 19466  9603  4813 11559 12927  8126 10522  1446
#> 2  13129     0 25491  8449 11894 24670 21309 15501 19012 11987
#> 3  19404 33675     0 18671 16446 19565 16158 12231 13538 19446
#> 4   9454  8232 18168     0  6929 15816 14800  8992 12503  9042
#> 5   4233 11909 16239  6149     0  9518 10098  5056  7841  3821
#> 6   7579 24609 19585 15946  8916     0  7604  8888  7902  7919
#> 7  11924 22042 15676 15262 10161  7477     0  7263  2989 11966
#> 8   7803 15674 12024  8893  5011  8698  8363     0  4351  7845
#> 9   9619 19259 15463 12479  7678  7696  3580  4381     0  9661
#> 10  1647 12015 20111  8626  3836  8413 12000  8683 10123     0

# a TSP tour based on durations:
tour <- 
  distA$durations |>
  as.dist() |> 
  TSP::TSP() |>
  TSP::solve_TSP()

tour
#> object of class 'TOUR' 
#> result of method 'arbitrary_insertion+two_opt' for 10 cities
#> tour length: 145

# which is basiaclly a named vector of point indeces
str(tour)
#>  'TOUR' Named int [1:10] 1 6 9 7 3 8 4 2 5 10
#>  - attr(*, "method")= chr "arbitrary_insertion+two_opt"
#>  - attr(*, "tour_length")= num 145
#>  - attr(*, "names")= chr [1:10] "1" "6" "9" "7" ...

# and can be converted to numeric or directly used for subsetting location frame
apotheke.df[tour,]
#>            id      lon      lat
#> 1   440338666 13.43853 52.47728
#> 6  3666540067 13.33844 52.49479
#> 9  1298511709 13.38544 52.53797
#> 7  4234162689 13.35448 52.54872
#> 3   977657079 13.49039 52.61719
#> 8  1492336686 13.43178 52.53216
#> 4  3770254015 13.51974 52.49737
#> 2   538057637 13.57874 52.45823
#> 5   364363337 13.45582 52.50112
#> 10  268726812 13.44518 52.47895

# and get a (single) route that passes all points in a sequence defined by `tour`
tour_sf <- osrmRoute(loc = apotheke.df[tour, c("lon", "lat")])
tour_sf
#> Simple feature collection with 1 feature and 4 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 13.33849 ymin: 52.45085 xmax: 13.58084 ymax: 52.62253
#> Geodetic CRS:  WGS 84
#>      src dst duration distance                       geometry
#> 1_10   1  10  143.905  84.1018 LINESTRING (13.43857 52.477...


leaflet(tour_sf) |> 
  addTiles()  |>
  addPolylines(label = ~ sprintf("%d min - %d km", round(duration), round(distance)),
               labelOptions = labelOptions(noHide = TRUE)) |>
  addMarkers(data = apotheke.df[tour,c("lon", "lat")])

If you need individual segments, it would still make sense to request only routes that are part of TSP results, perhaps something like this:

library(dplyr)
library(tidyr)

tour_segments_sf <- 
  apotheke.df[tour, ] |>
  # for destination coordinates, shift src coordinates by one row -- lead()
  mutate(across(c(lon,lat), lead, .names = "{.col}_dst")) |>
  na.omit() |>
  # osrmRoute only works with a single source, so we'll switch to rowwise grouping
  # (osrmRoute is called multiple times, once per each row)
  rowwise() |>
  mutate(route = osrmRoute(src = pick(lon, lat),
                           dst = pick(lon_dst, lat_dst))
  ) |>
  ungroup() |>
  unnest(route) |>
  st_as_sf()

tour_segments_sf |> select(id, duration, distance)
#> Simple feature collection with 9 features and 3 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 13.33849 ymin: 52.45085 xmax: 13.57897 ymax: 52.62253
#> Geodetic CRS:  WGS 84
#> # A tibble: 9 × 4
#>           id duration distance                                          geometry
#>        <dbl>    <dbl>    <dbl>                                  <LINESTRING [°]>
#> 1  440338666    15.0     11.6  (13.43857 52.47726, 13.44109 52.47274, 13.44235 …
#> 2 3666540067    14.6      7.90 (13.33849 52.4947, 13.34093 52.49981, 13.34624 5…
#> 3 1298511709     8.30     3.58 (13.38536 52.53791, 13.38357 52.5388, 13.38234 5…
#> 4 4234162689    24.2     15.7  (13.35467 52.54884, 13.35226 52.5504, 13.36443 5…
#> 5  977657079    20.9     12.2  (13.49051 52.61729, 13.49096 52.6167, 13.48765 5…
#> 6 1492336686    15.7      8.89 (13.43168 52.5322, 13.43243 52.53284, 13.4345 52…
#> 7 3770254015    15.2      8.23 (13.5196 52.4976, 13.52036 52.49818, 13.51961 52…
#> 8  538057637    20.0     11.9  (13.57897 52.45816, 13.57624 52.45456, 13.57197 …
#> 9  364363337     8.24     3.82 (13.45592 52.5011, 13.45536 52.5003, 13.46432 52…

leaflet(tour_segments_sf) |> 
  addTiles()  |>
  addPolylines(label = ~ sprintf("%d min - %d km", round(duration), round(distance)),
               labelOptions = labelOptions(noHide = TRUE)) |>
  addMarkers(data = apotheke.df[tour,c("lon", "lat")])

If it would not work for you, please add more details regarding your constraints and scope of your work.

Depends on the actual problem and topic, but for thesis, I would have assumed that http://project-osrm.org/ or any other remotely hosted routing service is not an option and you'd need to host your own osrm for experiments and/or work with OpenStreetMap data and network analysis tools (i.e. igraph, sfnetworks).