如何使用 data.table 有效地计算坐标对之间的距离：=

Question

我想找到最有效（最快）的方法来计算经纬度坐标对之间的距离。

使用sapply 和spDistsN1{sp} 提出了一个不太有效的解决方案(here)。我相信如果有人在data.table 中使用spDistsN1{sp} 和:= 运算符，这可以做得更快，但我无法做到这一点。有什么建议吗？

这是一个可重现的示例：

# load libraries
  library(data.table)
  library(dplyr)
  library(sp)
  library(rgeos)
  library(UScensus2000tract)

# load data and create an Origin-Destination matrix
  data("oregon.tract")

# get centroids as a data.frame
  centroids <- as.data.frame(gCentroid(oregon.tract,byid=TRUE))

# Convert row names into first column
  setDT(centroids, keep.rownames = TRUE)[]

# create Origin-destination matrix
  orig <- centroids[1:754, ]
  dest <- centroids[2:755, ]
  odmatrix <- bind_cols(orig,dest)
  colnames(odmatrix) <- c("origi_id", "long_orig", "lat_orig", "dest_id", "long_dest", "lat_dest")

我使用data.table 的尝试失败了

odmatrix[ , dist_km := spDistsN1(as.matrix(long_orig, lat_orig), as.matrix(long_dest, lat_dest), longlat=T)]

这是一个可行的解决方案（但可能效率较低）

odmatrix$dist_km <- sapply(1:nrow(odmatrix),function(i)
  spDistsN1(as.matrix(odmatrix[i,2:3]),as.matrix(odmatrix[i,5:6]),longlat=T))

head(odmatrix)

>   origi_id long_orig lat_orig  dest_id long_dest lat_dest dist_km
>      (chr)     (dbl)    (dbl)    (chr)     (dbl)    (dbl)   (dbl)
> 1 oregon_0   -123.51   45.982 oregon_1   -123.67   46.113 19.0909
> 2 oregon_1   -123.67   46.113 oregon_2   -123.95   46.179 22.1689
> 3 oregon_2   -123.95   46.179 oregon_3   -123.79   46.187 11.9014
> 4 oregon_3   -123.79   46.187 oregon_4   -123.83   46.181  3.2123
> 5 oregon_4   -123.83   46.181 oregon_5   -123.85   46.182  1.4054
> 6 oregon_5   -123.85   46.182 oregon_6   -123.18   46.066 53.0709

Answer 1

我编写了自己的 geosphere::distHaversine 版本，以便它更自然地适合 data.table := 调用，它可能在这里有用

dt.haversine <- function(lat_from, lon_from, lat_to, lon_to, r = 6378137){
  radians <- pi/180
  lat_to <- lat_to * radians
  lat_from <- lat_from * radians
  lon_to <- lon_to * radians
  lon_from <- lon_from * radians
  dLat <- (lat_to - lat_from)
  dLon <- (lon_to - lon_from)
  a <- (sin(dLat/2)^2) + (cos(lat_from) * cos(lat_to)) * (sin(dLon/2)^2)
  return(2 * atan2(sqrt(a), sqrt(1 - a)) * r)
}

2019 年 7 月 18 日更新

也可以通过 Rcpp 编写 C++ 版本。

#include <Rcpp.h>
using namespace Rcpp;

double inverseHaversine(double d){
  return 2 * atan2(sqrt(d), sqrt(1 - d)) * 6378137.0;
}

double distanceHaversine(double latf, double lonf, double latt, double lont,
                         double tolerance){
  double d;
  double dlat = latt - latf;
  double dlon =  lont - lonf;
  
  d = (sin(dlat * 0.5) * sin(dlat * 0.5)) + (cos(latf) * cos(latt)) * (sin(dlon * 0.5) * sin(dlon * 0.5));
  if(d > 1 && d <= tolerance){
    d = 1;
  }
  return inverseHaversine(d);
}

double toRadians(double deg){
  return deg * 0.01745329251;  // PI / 180;
}

// [[Rcpp::export]]
Rcpp::NumericVector rcpp_distance_haversine(Rcpp::NumericVector latFrom, Rcpp::NumericVector lonFrom, 
                        Rcpp::NumericVector latTo, Rcpp::NumericVector lonTo,
                        double tolerance) {
  
  int n = latFrom.size();
  NumericVector distance(n);
  
  double latf;
  double latt;
  double lonf;
  double lont;
  double dist = 0;
  
  for(int i = 0; i < n; i++){
    
    latf = toRadians(latFrom[i]);
    lonf = toRadians(lonFrom[i]);
    latt = toRadians(latTo[i]);
    lont = toRadians(lonTo[i]);
    dist = distanceHaversine(latf, lonf, latt, lont, tolerance);
    
    distance[i] = dist;
  }
  return distance;
}

将此文件保存在某处并使用Rcpp::sourceCpp("distance_calcs.cpp") 将函数加载到您的 R 会话中。

以下是它们与原始 geosphere::distHaversine 和 geosphere::distGeo 相比的一些基准

我将对象设置为 85k 行只是为了更有意义

dt <- rbindlist(list(odmatrix, odmatrix, odmatrix, odmatrix, odmatrix, odmatrix))
dt <- rbindlist(list(dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt, dt))

dt1 <- copy(dt); dt2 <- copy(dt); dt3 <- copy(dt); dt4 <- copy(dt)


library(microbenchmark)

microbenchmark(
  
  rcpp = {
    dt4[, dist := rcpp_distance_haversine(lat_orig, long_orig, lat_dest, long_dest, tolerance = 10000000000.0)]
  },
  
  dtHaversine = {
    dt1[, dist := dt.haversine(lat_orig, long_orig, lat_dest, long_dest)]
  }   ,
  
  haversine = {
    dt2[ , dist := distHaversine(matrix(c(long_orig, lat_orig), ncol = 2), 
                                 matrix(c(long_dest, lat_dest), ncol = 2))]
  },
  
  geo = {
    dt3[ , dist := distGeo(matrix(c(long_orig, lat_orig), ncol = 2), 
                           matrix(c(long_dest, lat_dest), ncol = 2))]
  },
  times = 5
)

# Unit: milliseconds
#       expr       min        lq      mean    median        uq        max neval
#        rcpp  5.622847  5.683959  6.208954  5.925277  6.036025   7.776664     5
# dtHaversine  9.024500 12.413380 12.335681 12.992920 13.590566  13.657037     5
#   haversine 30.911136 33.628153 52.503700 36.038927 40.791089 121.149197     5
#         geo 83.646104 83.971163 88.694377 89.548176 90.569327  95.737117     5

当然，由于两种不同技术（geo 和 hasrsine）计算距离的方式，结果会略有不同。

Answer 2

感谢@chinsoon12 的评论，我找到了一个结合distGeo{geosphere} 和data.table 的快速解决方案。在我的笔记本电脑中，快速解决方案比替代方案快 120 倍。

让我们把数据集做大来比较速度表现。

# Multiplicate data observations by 1000 
  odmatrix <- odmatrix[rep(seq_len(nrow(odmatrix)), 1000), ]

慢解

system.time(
           odmatrix$dist_km <- sapply(1:nrow(odmatrix),function(i)
             spDistsN1(as.matrix(odmatrix[i,2:3]),as.matrix(odmatrix[i,5:6]),longlat=T)) 
            )

 >   user  system elapsed 
 >   222.17    0.08  222.84 

快速解决方案

# load library
  library(geosphere)

# convert the data.frame to a data.table
  setDT(odmatrix)

system.time(
            odmatrix[ , dist_km2 := distGeo(matrix(c(long_orig, lat_orig), ncol = 2), 
                                            matrix(c(long_dest, lat_dest), ncol = 2))/1000]
           )

>   user  system elapsed 
>   1.76    0.03    1.79