Calculate, decode and plot routes on map using leaflet and R

Question

I have raw data which consists of lattitude and longitude of places The sample data is as follows:

EDIT (dput):

structure(list(Lat = c(-33.9409444, -33.9335713, -33.9333906, 
-33.9297826), Lon = c(18.5001774, 18.5033218, 18.518719, 18.5209372
)), .Names = c("Lat", "Lon"), row.names = c(NA, 4L), class = "data.frame")

I want to plot routes on the map using this data. This is my R code:

library(RODBC)
library(leaflet)

ui <- fluidPage(
  titlePanel("South Africa & Lesotho"),
  mainPanel(
    leafletOutput("mymap")
  )
)

server <- function(input, output, session) {
  dbhandle <- odbcDriverConnect('driver={SQL Server};server=localhost\\SQLEXpress;database=OSM;trusted_connection=true')
  res <- sqlQuery(dbhandle, 'select Lat, Lon from OSM2 where Street is not null')
  output$mymap <- renderLeaflet({
    leaflet(res) %>%
      addTiles() %>%
      addPolylines(lat = ~Lat, lng = ~Lon)
  }) 
}

shinyApp(ui, server)

However, all I get is this:

How can I use leaflet and R to plot the routes using the raw data (lat, long)?

What routes do you want to plot? `(lat, lon)` is only a position. `addPolylines` just connects one to another and that is what you see. — Christoph, Dec 09 '16 at 13:41
@Christoph I want to plot routes between the points given by the lat and lon — Munashe, Dec 09 '16 at 13:44
Please use `dput` when sharing your data. It is particularly difficult to reproduce a data frame from raw text when there are spaces in some of your text elements, such as in the `Street` column. — nrussell, Dec 09 '16 at 13:44
If you connect any point to any other point (with a real calculated route from e.g. osrm) your will see even less - everything will be blue. What is your goal? — Christoph, Dec 09 '16 at 13:58
@Christoph I want to plot the routes in a country (in this case South Africa). I extracted the raw data (lat, long etc) from osm and dumped it in an SQL Server database. So I now want to use this data to plot the routes in the country — Munashe, Dec 09 '16 at 14:02
Please `dput(data)` (or even better a small subset `dput(data[1:4, ]))` as already said... — Christoph, Dec 09 '16 at 14:09
@Munashe If the my answer is your solution, please mark the answer as solved as I put quite some effort into it;-) — Christoph, Dec 11 '16 at 12:55

Christoph · Accepted Answer · 2018-04-26T10:55:16.780

What you have to do:

Import the points
Calculate all routes between the points (I use OSRM)
Extract the route geometry from the routes (Appreciate the reference and have a look there for the speed updates!). Thanks to @SymbolixAU: You can also use googleway::decode_pl() or gepaf::decodePolyline()
Display everything on a map (I use leaflet)

My approach is not optimized for anything, but it should do the job... (It is script in RStudio, therefore the print() statements after leaflet.)

library(leaflet)
library(stringr)
library(bitops)

df <- structure(list(
  lat = c(-33.9409444, -33.9335713, -33.9333906, -33.9297826), 
  lng = c(18.5001774, 18.5033218, 18.518719, 18.5209372)),
  .Names = c("lat", "lng"), 
  row.names = c(NA, 4L), class = "data.frame")
nn <- nrow(df)

# Functions
# =========
viaroute <- function(lat1, lng1, lat2, lng2) {
  R.utils::evalWithTimeout({
    repeat {
      res <- try(
        route <- rjson::fromJSON(
          file = paste("http://router.project-osrm.org/route/v1/driving/",
                       lng1, ",", lat1, ";", lng2, ",", lat2,
                       "?overview=full", sep = "", NULL)))
      if (class(res) != "try-error") {
        if (!is.null(res)) {
          break
        }
      }
    }
  }, timeout = 1, onTimeout = "warning")
  return(res)
}

decode_geom <- function(encoded) {
  scale <- 1e-5
  len = str_length(encoded)
  encoded <- strsplit(encoded, NULL)[[1]]
  index = 1
  N <- 100000
  df.index <- 1
  array = matrix(nrow = N, ncol = 2)
  lat <- dlat <- lng <- dlnt <- b <- shift <- result <- 0

  while (index <= len) {
    # if (index == 80) browser()
    shift <- result <- 0
    repeat {
      b = as.integer(charToRaw(encoded[index])) - 63
      index <- index + 1
      result = bitOr(result, bitShiftL(bitAnd(b, 0x1f), shift))
      shift = shift + 5
      if (b < 0x20) break
    }
    dlat = ifelse(bitAnd(result, 1),
                  -(result - (bitShiftR(result, 1))),
                  bitShiftR(result, 1))
    lat = lat + dlat;

    shift <- result <- b <- 0
    repeat {
      b = as.integer(charToRaw(encoded[index])) - 63
      index <- index + 1
      result = bitOr(result, bitShiftL(bitAnd(b, 0x1f), shift))
      shift = shift + 5
      if (b < 0x20) break
    }
    dlng = ifelse(bitAnd(result, 1),
                  -(result - (bitShiftR(result, 1))),
                  bitShiftR(result, 1))
    lng = lng + dlng

    array[df.index,] <- c(lat = lat * scale, lng = lng * scale)
    df.index <- df.index + 1
  }

  geometry <- data.frame(array[1:df.index - 1,])
  names(geometry) <- c("lat", "lng")
  return(geometry)
}

map <- function() {
  m <- leaflet() %>%
    addTiles(group = "OSM") %>%
    addProviderTiles("Stamen.TonerLite") %>%
    addLayersControl(
      baseGroups = c("OSM", "Stamen.TonerLite")
    )
  return(m)
}

map_route <- function(df, my_list) {
  m <- map()
  m <- addCircleMarkers(map = m,
                        lat = df$lat,
                        lng = df$lng,
                        color = "blue",
                        stroke = FALSE,
                        radius = 6,
                        fillOpacity = 0.8) %>%
    addLayersControl(baseGroups = c("OSM", "Stamen.TonerLite")) %>%
    {
      for (i in 1:length(my_list)) {
        . <- addPolylines(., lat = my_list[[i]]$lat, lng = my_list[[i]]$lng, color = "red", weight = 4)
      }
      return(.)
    }
  return(m)
}

# Main
# ======
m <- map()
m <- m %>% addCircleMarkers(lat = df$lat,
                       lng = df$lng,
                       color = "red",
                       stroke = FALSE,
                       radius = 10,
                       fillOpacity = 0.8)
print(m)

my_list <- list()
r <- 1
for (i in 1:(nn-1)) {
  for (j in ((i+1):nn)) {
    my_route <- viaroute(df$lat[i], df$lng[i],df$lat[j], df$lng[j])
    geom <- decode_geom(my_route$routes[[1]]$geometry)
    my_list[[r]] <- geom
    r <- r + 1
  }
}

print(map_route(df, my_list))

Result:

In the end, you have to put all that in your shiny server...
I hope that helps!

if your `decode_geom`, is doing what I think it's doing (decoding polylines?) there is also `googleway::decode_pl()` and also `gepaf::decodePolyline()` — SymbolixAU, Dec 20 '16 at 21:19
and the benchmarks from [this answer](http://stackoverflow.com/a/32477722/5977215) — SymbolixAU, Dec 21 '16 at 00:59
@SymbolixAU thanks for pointing me to the updates! I made an edit to the answer... — Christoph, Dec 21 '16 at 07:41
No problem; he's going to add `googleway` into the benchmarks at some point too :) — SymbolixAU, Dec 21 '16 at 07:56

score 3 · Answer 2 · answered Mar 30 '20 at 22:42

Another more efficient way to calculate routes between points is with the osrm package: Interface Between R and the OpenStreetMap-Based Routing Service OSRM. Look at this example:

library(osrm)
library(leaflet)

df = data.frame(com = c("A", "B", "C"),
                lon = c(31.043515, 31.029080, 31.002896),
                lat = c(-29.778562, -29.795506, -29.836168),
                time = as.POSIXct(c("2020-03-18 07:56:59","2020-03-18 12:28:58","2020-03-18 18:24:52")))


trips <- osrmTrip(df, returnclass="sf")
trip <- trips[[1]]$trip

leaflet(trip) %>% 
  addProviderTiles("Stamen.TonerLite", group = "OSM") %>% 
  addPolylines() %>%
  addCircleMarkers(lat = df$lat,
                   lng = df$lon,
                   popup = paste(df$com,"-",format(df$time,"%H:%M:%S")),
                   color = "red",
                   stroke = FALSE,
                   radius = 8,
                   fillOpacity = 0.8)

score 2 · Answer 3 · answered Feb 13 '19 at 17:53

2

For the viaroute function created.

Use "R.utils::withTimeout" instead of "R.utils::evalWithTimeout" because that is now defunct.

I hope this helps

answered Feb 13 '19 at 17:53

user7779697

41
4

score 1 · Answer 4 · answered Sep 03 '22 at 09:23

@Christoph's code is fantastic - although some of the functions no longer work as originally written, presumably due to breaking changes in R over time.

As @user7779697 points out above, the viaroute() function needs be updated to accommodate code changes to R.utils that saw deprecation of evalWithTimeout, replacing it with withTimeout as follows:

R.utils::withTimeout()

I also ran into issues with the map_route() function, corrected by removing the braces from the internal for loop.

I've pasted the full updated code below which works with R Version 4.2.1 - I take no credit for this excellent work, only to get it back up and running with version changes:

library(leaflet)
library(stringr)
library(bitops)

df <- structure(list(
  lat = c(-33.9409444, -33.9335713, -33.9333906, -33.9297826), 
  lng = c(18.5001774, 18.5033218, 18.518719, 18.5209372)),
  .Names = c("lat", "lng"), 
  row.names = c(NA, 4L), class = "data.frame")
nn <- nrow(df)

# Functions
# =========
viaroute <- function(lat1, lng1, lat2, lng2) {
  R.utils::withTimeout({
    repeat {
      res <- try(
        route <- rjson::fromJSON(
          file = paste("http://router.project-osrm.org/route/v1/driving/",
                       lng1, ",", lat1, ";", lng2, ",", lat2,
                       "?overview=full", sep = "", NULL)))
      if (class(res) != "try-error") {
        if (!is.null(res)) {
          break
        }
      }
    }
  }, timeout = 1, onTimeout = "warning")
  return(res)
}

decode_geom <- function(encoded) {
  scale <- 1e-5
  len = str_length(encoded)
  encoded <- strsplit(encoded, NULL)[[1]]
  index = 1
  N <- 100000
  df.index <- 1
  array = matrix(nrow = N, ncol = 2)
  lat <- dlat <- lng <- dlnt <- b <- shift <- result <- 0
  
  while (index <= len) {
    # if (index == 80) browser()
    shift <- result <- 0
    repeat {
      b = as.integer(charToRaw(encoded[index])) - 63
      index <- index + 1
      result = bitOr(result, bitShiftL(bitAnd(b, 0x1f), shift))
      shift = shift + 5
      if (b < 0x20) break
    }
    dlat = ifelse(bitAnd(result, 1),
                  -(result - (bitShiftR(result, 1))),
                  bitShiftR(result, 1))
    lat = lat + dlat;
    
    shift <- result <- b <- 0
    repeat {
      b = as.integer(charToRaw(encoded[index])) - 63
      index <- index + 1
      result = bitOr(result, bitShiftL(bitAnd(b, 0x1f), shift))
      shift = shift + 5
      if (b < 0x20) break
    }
    dlng = ifelse(bitAnd(result, 1),
                  -(result - (bitShiftR(result, 1))),
                  bitShiftR(result, 1))
    lng = lng + dlng
    
    array[df.index,] <- c(lat = lat * scale, lng = lng * scale)
    df.index <- df.index + 1
  }
  
  geometry <- data.frame(array[1:df.index - 1,])
  names(geometry) <- c("lat", "lng")
  return(geometry)
}

map <- function() {
  m <- leaflet() %>%
    addTiles(group = "OSM") %>%
    addProviderTiles("Stamen.TonerLite") %>%
    addLayersControl(
      baseGroups = c("OSM", "Stamen.TonerLite")
    )
  return(m)
}

map_route <- function(df, my_list) {
  m <- map()
  m <- addCircleMarkers(map = m,
                        lat = df$lat,
                        lng = df$lng,
                        color = "blue",
                        stroke = FALSE,
                        radius = 6,
                        fillOpacity = 0.8) %>%
    addLayersControl(baseGroups = c("OSM", "Stamen.TonerLite")) 

      for (i in 1:length(my_list)) {
        m <- addPolylines(map = m, lat = my_list[[i]]$lat, lng = my_list[[i]]$lng, color = "red", weight = 4)
      }
      return(m)
}

# Main
# ======
m <- map()
m <- m %>% addCircleMarkers(lat = df$lat,
                            lng = df$lng,
                            color = "red",
                            stroke = FALSE,
                            radius = 10,
                            fillOpacity = 0.8)
print(m)


my_list <- list()
r <- 1
for (i in 1:(nn-1)) {
  for (j in ((i+1):nn)) {
    my_route <- viaroute(df$lat[i], df$lng[i],df$lat[j], df$lng[j])
    geom <- decode_geom(my_route$routes[[1]]$geometry)
    my_list[[r]] <- geom
    r <- r + 1
  }
}

print(map_route(df, my_list))

Calculate, decode and plot routes on map using leaflet and R

4 Answers4

Linked