使用leaflet和R计算、解码和绘制地图上的路线。

6

我将需要翻译有关IT技术的内容。以下是一些包含地点纬度和经度信息的原始数据样例:

编辑(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")

我希望使用这些数据在地图上绘制路线。这是我的 R 代码:

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)

然而,我得到的只有这个: enter image description here 如何使用leaflet和R绘制使用原始数据(lat,long)的路线?

1
你想绘制哪些路线?(lat, lon)只是一个位置。addPolylines只是将它们连接起来,这就是你看到的。 - Christoph
1
@Christoph 我想要在给定的纬度和经度点之间绘制路线。 - Munashe
1
请在分享数据时使用 dput。当您的文本元素中存在空格时,例如在“Street”列中,从原始文本中重现数据框架特别困难。 - nrussell
1
如果您将任何一个点与另一个点连接起来(例如使用osrm计算出真实的路线),您会看到更少的东西 - 一切都会变成蓝色。您的目标是什么? - Christoph
2
请执行 dput(data) (或者更好的是一个小子集 dput(data[1:4, ])),如前所述... - Christoph
显示剩余3条评论
4个回答

9

您需要做的:

  • 导入点
  • 计算所有点之间的路线(我使用OSRM
  • 从路线中提取路线几何形状(感谢参考并查看那里的速度更新!)。感谢@SymbolixAU:您还可以使用googleway :: decode_pl()gepaf :: decodePolyline()
  • 在地图上显示所有内容(我使用leaflet

我的方法没有针对任何优化,但应该能胜任...(这是RStudio中的脚本,因此在leaflet后使用print()语句。)

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

结果:

带有路线的点

最后,您需要将所有内容放入您的闪亮服务器中...
希望这能帮到您!


如果你的 decode_geom 函数正在执行我认为它在执行的操作(解码折线?),那么还有 googleway::decode_pl()gepaf::decodePolyline() 可以使用。 - SymbolixAU
1
返回以下答案中的基准测试结果:https://dev59.com/v47ea4cB1Zd3GeqPHu1Q#32477722 - SymbolixAU
@SymbolixAU谢谢您指出更新!我对答案进行了编辑... - Christoph
没问题,他会在某个时候将 googleway 加入基准测试中 :) - SymbolixAU

3

另一种更高效的计算点之间路径的方法是使用osrm包:它是R和基于OpenStreetMap的路由服务OSRM之间的接口。看一个例子:

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


2

对于创建的"viaroute"函数。

请使用"R.utils::withTimeout"而不是"R.utils::evalWithTimeout",因为后者现已废弃。

希望这可以帮到你。

注:Original Answer翻译成"最初的回答"。


1

@Christoph的代码非常棒 - 尽管一些函数不再像最初编写的那样工作,这可能是由于R随着时间的推移而出现了破坏性变化。

正如@user7779697在上面指出的那样,viaroute()函数需要更新以适应对R.utils代码更改的调整,该代码看到evalWithTimeout被弃用,并将其替换为withTimeout,如下所示:

R.utils::withTimeout()

我也遇到了与 map_route() 函数有关的问题,通过从内部 for 循环中删除括号进行更正。

我已经粘贴下面完整的更新代码,它可以在 R 版本 4.2.1 中运行 - 我对这个出色的工作没有任何功劳,只是为了让它适应版本更改而使其重新运行:

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

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接