如何解码来自 OSRM 的编码折线并绘制路线几何图形？
How to decode encoded polylines from OSRM and plotting route geometry?
我正在使用 OSRM（OpenStreetMap Routing Machine）实例来评估不同点的距离和时间。使用 API，我可以检索我想要和需要的信息，尤其是作为折线的真实路线。
直到今天，我已经绘制了起点和终点之间的直线。
segments(
lon_patient,lat_patient,lon_lieu,lat_lieu,col = transp_time,lwd = 3
)
现在我想绘制多段线。但它是编码的（https://github.com/Project-OSRM/osrm-backend/wiki/Server-api#response-2）。怎么画呢？
谢谢！
一种（快速）方法是从 mapbox github 存储库下载 polyline.js 文件，然后使用 V8 包为您完成艰苦的工作：
library(V8)
ctx <- new_context()
ctx$source("polyline.js")
ctx$call("polyline.decode", "_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## [,1] [,2]
## [1,] 38.500 -120.200
## [2,] 40.700 -120.950
## [3,] 43.252 -126.453
它 returns 一个 lat/lon 对矩阵，你应该可以使用它。
纯 R/Rcpp 答案在长 运行 中会更好。
更新
有一个！这来自：https://gist.github.com/diegovalle/916889（我添加了 require 并组合了一些冗长的 0 作业）：
DecodeLineR <- function(encoded) {
require(bitops)
require(stringr)
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) {
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 * 1e-05, lng = lng * 1e-5)
df.index <- df.index + 1
}
ret <- data.frame(array[1:df.index - 1,])
names(ret) <- c("lat", "lng")
return(ret)
}
DecodeLineR("_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## lat lng
## 1 38.500 -120.200
## 2 40.700 -120.950
## 3 43.252 -126.453
这会让你得到一个数据框和一个矩阵。并且是纯 R。不确定哪个会更快（如果需要速度）。
更新 #2
这里还有另一个纯 R 实现：http://s4rdd.blogspot.com/2012/12/google-maps-api-decoding-polylines-for.html 它比上面的要快得多（参见下面的基准测试）。
decodeLine <- function(encoded){
require(bitops)
vlen <- nchar(encoded)
vindex <- 0
varray <- NULL
vlat <- 0
vlng <- 0
while(vindex < vlen){
vb <- NULL
vshift <- 0
vresult <- 0
repeat{
if(vindex + 1 <= vlen){
vindex <- vindex + 1
vb <- as.integer(charToRaw(substr(encoded, vindex, vindex))) - 63
}
vresult <- bitOr(vresult, bitShiftL(bitAnd(vb, 31), vshift))
vshift <- vshift + 5
if(vb < 32) break
}
dlat <- ifelse(
bitAnd(vresult, 1)
, -(bitShiftR(vresult, 1)+1)
, bitShiftR(vresult, 1)
)
vlat <- vlat + dlat
vshift <- 0
vresult <- 0
repeat{
if(vindex + 1 <= vlen) {
vindex <- vindex+1
vb <- as.integer(charToRaw(substr(encoded, vindex, vindex))) - 63
}
vresult <- bitOr(vresult, bitShiftL(bitAnd(vb, 31), vshift))
vshift <- vshift + 5
if(vb < 32) break
}
dlng <- ifelse(
bitAnd(vresult, 1)
, -(bitShiftR(vresult, 1)+1)
, bitShiftR(vresult, 1)
)
vlng <- vlng + dlng
varray <- rbind(varray, c(vlat * 1e-5, vlng * 1e-5))
}
coords <- data.frame(varray)
names(coords) <- c("lat", "lon")
coords
}
这是 Rcpp/C++11 版本，由 https://mapzen.com/documentation/mobility/decoding/ 提供：
#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
DataFrame decode_polyline(const std::string& encoded) {
size_t i = 0; // what byte are we looking at
constexpr double kPolylinePrecision = 1E6;
constexpr double kInvPolylinePrecision = 1.0 / kPolylinePrecision;
auto deserialize = [&encoded, &i](const int previous) {
int byte, shift = 0, result = 0;
do {
byte = static_cast<int>(encoded[i++]) - 63;
result |= (byte & 0x1f) << shift;
shift += 5;
} while (byte >= 0x20);
return previous + (result & 1 ? ~(result >> 1) : (result >> 1));
};
std::vector<double> lonv, latv;
int last_lon = 0, last_lat = 0;
while (i < encoded.length()) {
int lat = deserialize(last_lat);
int lon = deserialize(last_lon);
latv.emplace_back(static_cast<float>(static_cast<double>(lat) * kInvPolylinePrecision));
lonv.emplace_back(static_cast<float>(static_cast<double>(lon) * kInvPolylinePrecision));
last_lon = lon;
last_lat = lat;
}
return DataFrame::create(_["lon"] = lonv, _["lat"] = latv);
}
将其保存到 polyline.cpp，然后：
Rcpp::sourceCpp("polyline.cpp")
那么您可以：
decode_polyline("_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## lon lat
## 1 -120.200 38.500
## 2 -120.950 40.700
#3 3 -126.453 43.252
基准
我将这两个 R 函数采购到全局环境中，并为 javascript 和 C++ 实现执行了 js 和 C++ 等价物。
无论我使用什么微基准参数，DecodeLineR 的最大值都相当 "out there"。 decodeLine() 纯 R 版本的性能似乎足以保证不会产生 V8 或 Rcpp/C++11 依赖性，但是 YMMV.
最终更新（MOAR 基准）
我将 googleway::decode_pl() 函数合并到新的基准测试中，并使用了更长的折线。基准代码在下面，新情节在下面。
library(microbenchmark)
library(Rcpp)
library(inline)
library(V8)
library(googleway)
library(ggplot2)
sourceCpp("polyline.cpp")
ctx <- v8()
ctx$source("polyline.js")
source("DecodeLineR.R")
source("decodeline.R")
line_str <- "{ae{HntiQtCcDzG_I|^uc@rFgHhC{CxAiA~AaA~BkAvB}A|F_G|AgBbBkCtAwCd@sA|BoIVw@Pc@|@gBt@}@|@y@lCwBvA_B`@k@~@aBt@iBlAaE~@oEp@sDX{BP_BJaDAcEIeCe@gHo@yMUaEk@uDm@iD]mCAwBNsDXyDL}@nByIZyCt@cLr@gNB_ABoEAkFJmDTkBVeAZ_Af@gAnDwF|@gBbAoChHgUPWlAT`@B|@GbE_@dAW`Cu@vBe@tDs@xD{@`Bg@bBq@hBaAtB}@dCi@bF}@jBg@pBeAj@SNE\C^@\DbAZ`Ah@~C`A\H|ALzAFLA^Gl@UdBgAjBaBZSh@Qz@MjD_@`FoAtCa@j@Ez@DxE|@xF\nBP~@TxHvBf@Tb@\pBvC\^`@XxAf@fBT|BDfAIr@MfBe@rBa@rBMvBYxBg@xA_@^Ir@@NF|@l@nBfAjAj@dBV`Bb@lBbAbB~ALPhC`FV`@n@z@^VNBX?LGZa@d@eAp@qAt@Sx@Cz@G\IZOhCcBb@c@T]jA_CrE_HfEiFz@}@p@k@|@o@`C{A`A{@rBwBx@oAbByCp@wArAoDLWxA}BhAcBjAqAlAiB~AaDr@sBp@{CD[TkC^}FZyD^oCx@gF`@qAh@kAz@yAtAgBpD_E|JoKdDuEjBcCfC{ExCqGdAgBlBuBrAyBpEkIpEsI\]^YbAg@|GaBzKeEfBe@lCW`AQr@U|A_AtAkAhDyCpAeA|Aq@`EeCrDgBvA{@tD}C`BmAzBm@t@QvAQxBOl@Q~Ai@~BsAlCcB"
microbenchmark(
googleway = decode_pl(line_str),
rcpp = decode(line_str),
js = ctx$call("polyline_decode", line_str),
DecodeLineR = DecodeLineR(line_str),
decodeLine = decodeLine(line_str),
control=list(warmup=50),
times=1000
) -> mb
mb
## Unit: microseconds
## expr min lq mean median uq max neval cld
## googleway 404.322 471.8475 817.8312 526.270 579.095 135564.973 1000 a
## rcpp 253.062 302.9550 363.9325 359.765 401.054 831.699 1000 a
## js 2793.677 3099.3390 3359.6190 3219.312 3427.636 15580.609 1000 b
## DecodeLineR 9366.714 9656.4140 12029.3991 10128.676 12570.216 181536.940 1000 c
## decodeLine 9907.869 10218.0465 11413.5732 10655.949 11491.305 150588.501 1000 c
update_geom_defaults("violin", list(fill="maroon"))
autoplot(mb) +
scale_y_log10(name="Time [microseconds]", label=scales::comma) +
hrbrmisc::theme_hrbrmstr(grid="X")
我正在使用 OSRM（OpenStreetMap Routing Machine）实例来评估不同点的距离和时间。使用 API，我可以检索我想要和需要的信息，尤其是作为折线的真实路线。
直到今天，我已经绘制了起点和终点之间的直线。
segments(
lon_patient,lat_patient,lon_lieu,lat_lieu,col = transp_time,lwd = 3
)
现在我想绘制多段线。但它是编码的（https://github.com/Project-OSRM/osrm-backend/wiki/Server-api#response-2）。怎么画呢？
谢谢！
一种（快速）方法是从 mapbox github 存储库下载 polyline.js 文件，然后使用 V8 包为您完成艰苦的工作：
library(V8)
ctx <- new_context()
ctx$source("polyline.js")
ctx$call("polyline.decode", "_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## [,1] [,2]
## [1,] 38.500 -120.200
## [2,] 40.700 -120.950
## [3,] 43.252 -126.453
它 returns 一个 lat/lon 对矩阵，你应该可以使用它。
纯 R/Rcpp 答案在长 运行 中会更好。
更新
有一个！这来自：https://gist.github.com/diegovalle/916889（我添加了 require 并组合了一些冗长的 0 作业）：
DecodeLineR <- function(encoded) {
require(bitops)
require(stringr)
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) {
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 * 1e-05, lng = lng * 1e-5)
df.index <- df.index + 1
}
ret <- data.frame(array[1:df.index - 1,])
names(ret) <- c("lat", "lng")
return(ret)
}
DecodeLineR("_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## lat lng
## 1 38.500 -120.200
## 2 40.700 -120.950
## 3 43.252 -126.453
这会让你得到一个数据框和一个矩阵。并且是纯 R。不确定哪个会更快（如果需要速度）。
更新 #2
这里还有另一个纯 R 实现：http://s4rdd.blogspot.com/2012/12/google-maps-api-decoding-polylines-for.html 它比上面的要快得多（参见下面的基准测试）。
decodeLine <- function(encoded){
require(bitops)
vlen <- nchar(encoded)
vindex <- 0
varray <- NULL
vlat <- 0
vlng <- 0
while(vindex < vlen){
vb <- NULL
vshift <- 0
vresult <- 0
repeat{
if(vindex + 1 <= vlen){
vindex <- vindex + 1
vb <- as.integer(charToRaw(substr(encoded, vindex, vindex))) - 63
}
vresult <- bitOr(vresult, bitShiftL(bitAnd(vb, 31), vshift))
vshift <- vshift + 5
if(vb < 32) break
}
dlat <- ifelse(
bitAnd(vresult, 1)
, -(bitShiftR(vresult, 1)+1)
, bitShiftR(vresult, 1)
)
vlat <- vlat + dlat
vshift <- 0
vresult <- 0
repeat{
if(vindex + 1 <= vlen) {
vindex <- vindex+1
vb <- as.integer(charToRaw(substr(encoded, vindex, vindex))) - 63
}
vresult <- bitOr(vresult, bitShiftL(bitAnd(vb, 31), vshift))
vshift <- vshift + 5
if(vb < 32) break
}
dlng <- ifelse(
bitAnd(vresult, 1)
, -(bitShiftR(vresult, 1)+1)
, bitShiftR(vresult, 1)
)
vlng <- vlng + dlng
varray <- rbind(varray, c(vlat * 1e-5, vlng * 1e-5))
}
coords <- data.frame(varray)
names(coords) <- c("lat", "lon")
coords
}
这是 Rcpp/C++11 版本，由 https://mapzen.com/documentation/mobility/decoding/ 提供：
#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
DataFrame decode_polyline(const std::string& encoded) {
size_t i = 0; // what byte are we looking at
constexpr double kPolylinePrecision = 1E6;
constexpr double kInvPolylinePrecision = 1.0 / kPolylinePrecision;
auto deserialize = [&encoded, &i](const int previous) {
int byte, shift = 0, result = 0;
do {
byte = static_cast<int>(encoded[i++]) - 63;
result |= (byte & 0x1f) << shift;
shift += 5;
} while (byte >= 0x20);
return previous + (result & 1 ? ~(result >> 1) : (result >> 1));
};
std::vector<double> lonv, latv;
int last_lon = 0, last_lat = 0;
while (i < encoded.length()) {
int lat = deserialize(last_lat);
int lon = deserialize(last_lon);
latv.emplace_back(static_cast<float>(static_cast<double>(lat) * kInvPolylinePrecision));
lonv.emplace_back(static_cast<float>(static_cast<double>(lon) * kInvPolylinePrecision));
last_lon = lon;
last_lat = lat;
}
return DataFrame::create(_["lon"] = lonv, _["lat"] = latv);
}
将其保存到 polyline.cpp，然后：
Rcpp::sourceCpp("polyline.cpp")
那么您可以：
decode_polyline("_p~iF~ps|U_ulLnnqC_mqNvxq`@")
## lon lat
## 1 -120.200 38.500
## 2 -120.950 40.700
#3 3 -126.453 43.252
基准
我将这两个 R 函数采购到全局环境中，并为 javascript 和 C++ 实现执行了 js 和 C++ 等价物。
无论我使用什么微基准参数，DecodeLineR 的最大值都相当 "out there"。 decodeLine() 纯 R 版本的性能似乎足以保证不会产生 V8 或 Rcpp/C++11 依赖性，但是 YMMV.
最终更新（MOAR 基准）
我将 googleway::decode_pl() 函数合并到新的基准测试中，并使用了更长的折线。基准代码在下面，新情节在下面。
library(microbenchmark)
library(Rcpp)
library(inline)
library(V8)
library(googleway)
library(ggplot2)
sourceCpp("polyline.cpp")
ctx <- v8()
ctx$source("polyline.js")
source("DecodeLineR.R")
source("decodeline.R")
line_str <- "{ae{HntiQtCcDzG_I|^uc@rFgHhC{CxAiA~AaA~BkAvB}A|F_G|AgBbBkCtAwCd@sA|BoIVw@Pc@|@gBt@}@|@y@lCwBvA_B`@k@~@aBt@iBlAaE~@oEp@sDX{BP_BJaDAcEIeCe@gHo@yMUaEk@uDm@iD]mCAwBNsDXyDL}@nByIZyCt@cLr@gNB_ABoEAkFJmDTkBVeAZ_Af@gAnDwF|@gBbAoChHgUPWlAT`@B|@GbE_@dAW`Cu@vBe@tDs@xD{@`Bg@bBq@hBaAtB}@dCi@bF}@jBg@pBeAj@SNE\C^@\DbAZ`Ah@~C`A\H|ALzAFLA^Gl@UdBgAjBaBZSh@Qz@MjD_@`FoAtCa@j@Ez@DxE|@xF\nBP~@TxHvBf@Tb@\pBvC\^`@XxAf@fBT|BDfAIr@MfBe@rBa@rBMvBYxBg@xA_@^Ir@@NF|@l@nBfAjAj@dBV`Bb@lBbAbB~ALPhC`FV`@n@z@^VNBX?LGZa@d@eAp@qAt@Sx@Cz@G\IZOhCcBb@c@T]jA_CrE_HfEiFz@}@p@k@|@o@`C{A`A{@rBwBx@oAbByCp@wArAoDLWxA}BhAcBjAqAlAiB~AaDr@sBp@{CD[TkC^}FZyD^oCx@gF`@qAh@kAz@yAtAgBpD_E|JoKdDuEjBcCfC{ExCqGdAgBlBuBrAyBpEkIpEsI\]^YbAg@|GaBzKeEfBe@lCW`AQr@U|A_AtAkAhDyCpAeA|Aq@`EeCrDgBvA{@tD}C`BmAzBm@t@QvAQxBOl@Q~Ai@~BsAlCcB"
microbenchmark(
googleway = decode_pl(line_str),
rcpp = decode(line_str),
js = ctx$call("polyline_decode", line_str),
DecodeLineR = DecodeLineR(line_str),
decodeLine = decodeLine(line_str),
control=list(warmup=50),
times=1000
) -> mb
mb
## Unit: microseconds
## expr min lq mean median uq max neval cld
## googleway 404.322 471.8475 817.8312 526.270 579.095 135564.973 1000 a
## rcpp 253.062 302.9550 363.9325 359.765 401.054 831.699 1000 a
## js 2793.677 3099.3390 3359.6190 3219.312 3427.636 15580.609 1000 b
## DecodeLineR 9366.714 9656.4140 12029.3991 10128.676 12570.216 181536.940 1000 c
## decodeLine 9907.869 10218.0465 11413.5732 10655.949 11491.305 150588.501 1000 c
update_geom_defaults("violin", list(fill="maroon"))
autoplot(mb) +
scale_y_log10(name="Time [microseconds]", label=scales::comma) +
hrbrmisc::theme_hrbrmstr(grid="X")