r - r - 按组分的地理距离 - 在每对行上应用一个函数

我想计算每个省的房子之间的平均地理距离。

假设我有以下数据。


df1 <- data.frame(province = c(1, 1, 1, 2, 2, 2),


 house = c(1, 2, 3, 4, 5, 6),


 lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 


 lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))



使用 geosphere 库,我可以找到两个房子之间的距离。 举个例子:


library(geosphere)


distm(c(df1$lon[1], df1$lat[1]), c(df1$lon[2], df1$lat[2]), fun = distHaversine)



#11429.1



如何计算省内各房屋间的距离并收集每个省的平均距离?

原始数据集每个省份有数百万个观测,所以性能也是一个问题。

时间:

解决方案:


lapply(split(df1, df1$province), function(df){


 df <- Expand.Grid(df[, c("lat","lon")], df[, c("lat","lon")])


 mean(distHaversine(df[, 1:2], df[, 3:4]))


})



其中 Expand.Grid() 是从 https://stackoverflow.com/a/30085602/3502164 获取的。

解释:

1.性能测试

在将 vectorised函数 distHaversine() 转换为 unvectorised distm() 时,我将避免使用 distm() 如果你查看你看到的源代码:


function (x, y, fun = distHaversine) 


{


 [...]


 for (i in 1:n) {


 dm[i, ] = fun(x[i, ], y)


 }


 return(dm)


}



distHaversine() 发送"整个对象"到C 时,distm() 将数据"行智能"发送到 distHaversine(),因这里在执行C 时强制 distHaversine() 执行相同的操作。 因此,不应该使用 distm() 。 在性能方面,我看到使用包装函数 distm() 带来更多危害,因为我看到。

2.解释"解决方案"中的代码:

a ) 分组:

要分析每个组的数据,请执行以下操作: province. Splitting可以通过以下方式完成: split(df1, df1$province)

b ) 分组"columns"

你想查找纬度/经度的所有唯一组合。 第一个猜测可能是 Expand.Grid(),但在多个列中不适用。 Flick先生在中负责这个 expand.grid 函数的data.frames 。

然后你就有了所有可能组合的data.frame(),只需要使用 mean(distHaversine(...))

我的初始想法是查看 distHaversine的源代码,并将它的复制到 proxy that的函数中。 就像这个( 请注意,lon 将是第一列):


library(geosphere)


library(dplyr)


library(proxy)



df1 <- data.frame(province = as.integer(c(1, 1, 1, 2, 2, 2)),


 house = as.integer(c(1, 2, 3, 4, 5, 6)),


 lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7), 


 lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2))



custom_haversine <- function(x, y) {


 toRad <- pi/180



 diff <- (y - x) * toRad


 dLon <- diff[1L]


 dLat <- diff[2L]



 a <- sin(dLat/2) ^ 2 + cos(x[2L] * toRad) * cos(y[2L] * toRad) * sin(dLon/2) ^ 2


 a <- min(a, 1)


 # return


 2 * atan2(sqrt(a), sqrt(1 - a)) * 6378137


}



pr_DB$set_entry(FUN=custom_haversine, names="haversine", loop=TRUE, distance=TRUE)



average_dist <- df1 %>%


 select(-house) %>%


 group_by(province) %>%


 group_map(~ data.frame(avg=mean(proxy::dist(.x[, c("lon","lat")], method="haversine"))))



但是,如果你期望每省数百行,proxy 可以能无法分配中间( 下三角形) 矩阵。 所以我把代码移植到 C++ 中,增加了多线程作为额外的奖励:

编辑:事实证明s2d助手远非最佳,这个版本现在使用这里给出的公式。


//[[Rcpp::plugins(cpp11)]]


//[[Rcpp::depends(RcppParallel)]]



#include <cstddef>//size_t


#include <math.h>//sin, cos, sqrt, atan2, pow


#include <vector>



#include <Rcpp.h>


#include <RcppParallel.h>



using namespace std;


using namespace Rcpp;


using namespace RcppParallel;



//single to double indices for lower triangular of matrices without diagonal


void s2d(const size_t id, const size_t nrow, size_t& i, size_t& j) {


 j = nrow - 2 - static_cast<size_t>(sqrt(-8 * id + 4 * nrow * (nrow - 1) - 7)/2 - 0.5);


 i = id + j + 1 - nrow * (nrow - 1)/2 + (nrow - j) * ((nrow - j) - 1)/2;


}



class HaversineCalculator : public Worker


{


public:


 HaversineCalculator(const NumericVector& lon,


 const NumericVector& lat,


 double& avg,


 const int n)


 : lon_(lon)


, lat_(lat)


, avg_(avg)


, n_(n)


, cos_lat_(lon.length())


 {


//terms for distance calculation


 for (size_t i = 0; i <cos_lat_.size(); i++) {


 cos_lat_[i] = cos(lat_[i] * 3.1415926535897/180);


 }


 }



 void operator()(size_t begin, size_t end) {


//for Kahan summation


 double sum = 0;


 double c = 0;



 double to_rad = 3.1415926535897/180;



 size_t i, j;


 for (size_t ind = begin; ind <end; ind++) {


 s2d(ind, lon_.length(), i, j);



//haversine distance


 double d_lon = (lon_[j] - lon_[i]) * to_rad;


 double d_lat = (lat_[j] - lat_[i]) * to_rad;


 double d_hav = pow(sin(d_lat/2), 2) + cos_lat_[i] * cos_lat_[j] * pow(sin(d_lon/2), 2);


 if (d_hav> 1) d_hav = 1;


 d_hav = 2 * atan2(sqrt(d_hav), sqrt(1 - d_hav)) * 6378137;



//the average part


 d_hav/= n_;



//Kahan sum step


 double y = d_hav - c;


 double t = sum + y;


 c = (t - sum) - y;


 sum = t;


 }



 mutex_.lock();


 avg_ += sum;


 mutex_.unlock();


 }



private:


 const RVector<double> lon_;


 const RVector<double> lat_;


 double& avg_;


 const int n_;


 tthread::mutex mutex_;


 vector<double> cos_lat_;


};



//[[Rcpp::export]]


double avg_haversine(const DataFrame& input, const int nthreads) {


 NumericVector lon = input["lon"];


 NumericVector lat = input["lat"];



 double avg = 0;


 int size = lon.length() * (lon.length() - 1)/2;


 HaversineCalculator hc(lon, lat, avg, size);



 int grain = size/nthreads/10;//play with this value for load balance


 parallelFor(0, size, hc, grain);



 return avg;


}



这个代码不会分配任何中间矩阵,它将简单地计算每一对的下三角形的距离。 请参见这里的,以了解求和部分。

如果将该代码保存为 haversine.cpp,则可以执行以下操作:


library(dplyr)


library(Rcpp)


library(RcppParallel)



sourceCpp("haversine.cpp")



df1 %>%


 group_by(province) %>%


 group_map(~ data.frame(avg=avg_haversine(.x, parallel::detectCores())))


# A tibble: 2 x 2


# Groups: province [2]


 province avg


 <int> <dbl>


1 1 15379.


2 2 793612.



下面是一个健全检查:


pr_DB$set_entry(FUN=geosphere::distHaversine, names="distHaversine", loop=TRUE, distance=TRUE)



df1 %>%


 select(-house) %>%


 group_by(province) %>%


 group_map(~ data.frame(avg=mean(proxy::dist(.x[, c("lon","lat")], method="distHaversine"))))



但要注意的是:


df <- data.frame(lon=runif(1e3, -90, 90), lat=runif(1e3, -90, 90))



system.time(proxy::dist(df, method="distHaversine"))


 user system elapsed 


 34.353 0.005 34.394



system.time(proxy::dist(df, method="haversine"))


 user system elapsed 


 0.789 0.020 0.809



system.time(avg_haversine(df, 4L))


 user system elapsed 


 0.054 0.000 0.014



df <- data.frame(lon=runif(1e5, -90, 90), lat=runif(1e5, -90, 90))



system.time(avg_haversine(df, 4L))


 user system elapsed 


 73.861 0.238 19.670



如果你有数百万行,你可能需要等待一段时间......

此外,我还应该提到,如果你通过创建的线程是不可能的,那么如果你开始计算,那么应该等到完成计算,或者重启 R/rstudio 。

关于复杂性

根据你的实际数据和计算机有多少内核,你可以能会很好地完成计算的计算。 这个问题有二次复杂性(每个省,可以这么说)。 这一行:


int size = lon.length() * (lon.length() - 1)/2;



表示必须执行的( haversine ) 距离计算的数量。 因这里,如果行数由 n的数量增加,计算的数量会随着 n^2/2的数量增加,大致提高。

如果你先计算每个数字,也不可能计算出平均值,并且你要在一台机器上运行一个或者多个机器,这样你就得花一段时间去寻找比多线程代码更快的东西。 否则你就不能解决这个问题。

...