Search code examples
rr-leaflet

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


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:

enter image description here

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


Solution

  • 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:

    Points with routes

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