CODEWITHSUNDEEP
rr-leaflet
Munashe
Munashe

Reputation: 101

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)?

Upvotes: 6

Views: 5390

Answers (4)

Nick
Nick

Reputation: 286

@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))

Upvotes: 1

Rafael D&#237;az
Rafael D&#237;az

Reputation: 2289

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)

enter image description here

Upvotes: 3

user7779697
user7779697

Reputation: 41

For the viaroute function created.

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

I hope this helps

Upvotes: 2

Christoph
Christoph

Reputation: 7063

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!

Upvotes: 9

Related Questions