按组划分的地理距离 - 在每对行上应用一个函数
Geographical distance by group - Applying a function on each pair of rows
我想计算每个省多个房屋之间的平均地理距离。
假设我有以下数据。
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
如何计算省内所有房屋之间的距离并收集每个省的平均距离？
原始数据集每个省都有数百万个观测值，因此性能也是一个问题。
我的 10 美分。您可以：
# subset the province
df1 <- df1[which(df1$province==1),]
# get all combinations
all <- combn(df1$house, 2, FUN = NULL, simplify = TRUE)
# run your function and get distances for all combinations
distances <- c()
for(col in 1:ncol(all)) {
a <- all[1, col]
b <- all[2, col]
dist <- distm(c(df1$lon[a], df1$lat[a]), c(df1$lon[b], df1$lat[b]), fun = distHaversine)
distances <- c(distances, dist)
}
# calculate mean:
mean(distances)
# [1] 15379.21
这为您提供了该省的平均值，您可以将其与其他方法的结果进行比较。例如评论中提到的sapply：
df1 <- df1[which(df1$province==1),]
mean(sapply(split(df1, df1$province), dist))
# [1] 1.349036
如您所见，它给出了不同的结果，因为 dist 函数可以计算不同类型的距离（如欧几里得）而不能计算半正弦或其他 "geodesic" 距离。包 geodist 似乎有选项可以让你更接近 sapply:
library(geodist)
library(magrittr)
# defining the data
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))
# defining the function
give_distance <- function(resultofsplit){
distances <- c()
for (i in 1:length(resultofsplit)){
sdf <- resultofsplit
sdf <- sdf[[i]]
sdf <- sdf[c("lon", "lat", "province", "house")]
sdf2 <- as.matrix(sdf)
sdf3 <- geodist(x=sdf2, measure="haversine")
sdf4 <- unique(as.vector(sdf3))
sdf4 <- sdf4[sdf4 != 0] # this is to remove the 0-distances
mean_dist <- mean(sdf4)
distances <- c(distances, mean_dist)
}
return(distances)
}
split(df1, df1$province) %>% give_distance()
#[1] 15379.21 793612.04
例如该函数将为您提供每个省的平均距离值。现在，我没有设法让 give_distance 与 sapply 一起工作，但这应该已经更有效率了。
参考此 ，您的问题的矢量化解决方案如下所示；
toCheck <- sapply(split(df1, df1$province), function(x){
combn(rownames(x), 2, simplify = FALSE)})
names(toCheck) <- sapply(toCheck, paste, collapse = " - ")
sapply(toCheck, function(x){
distm(df1[x[1],c("lon","lat")], df1[x[2],c("lon","lat")],
fun = distHaversine)
})
# 1 - 2 1 - 3 2 - 3 4 - 5 4 - 6 5 - 6
# 11429.10 22415.04 12293.48 634549.20 1188925.65 557361.28
如果每个省份的记录数相同，则此方法有效。如果不是这种情况，那么为 toCheck 分配适当名称的第二部分以及我们在最后如何使用它应该随着 toCheck 列表结构的变化而改变。它不关心数据集的顺序。
对于您的实际数据集，toCheck 将成为嵌套列表，因此您需要像下面这样调整函数；对于此解决方案，我还没有为 toCheck 命名。 （df2可以在答案末尾找到）。
df2 <- df2[order(df2$province),] #sorting may even improve performance
names(toCheck) <- paste("province", unique(df2$province))
toCheck <- sapply(split(df2, df2$province), function(x){
combn(rownames(x), 2, simplify = FALSE)})
sapply(toCheck, function(x){ sapply(x, function(y){
distm(df2[y[1],c("lon","lat")], df2[y[2],c("lon","lat")], fun = distHaversine)
})})
# $`province 1`
# [1] 11429.10 22415.04 1001964.84 12293.48 1013117.36 1024209.46
#
# $`province 2`
# [1] 634549.2 1188925.7 557361.3
#
# $`province 3`
# [1] 590083.2
#
# $`province 4`
# [1] 557361.28 547589.19 11163.92
您可以进一步获得每个省份的mean()。此外，如果需要，重命名嵌套列表的元素应该不难，这样您就可以知道每个距离对应于哪些房屋。
df2 <- data.frame(province = c(1, 1, 1, 2, 2, 2, 1, 3, 3, 4,4,4),
house = c(1, 2, 3, 4, 5, 6, 7, 10, 9, 8, 11, 12),
lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7, -85.6, -76.4, -75.4, -80.9, -85.7, -85.6),
lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2, 40.1, 39.3, 60.8, 53.3, 40.2, 40.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() 取自 。
解释：
1.性能
我会避免使用 distm()，因为它将 vectorised 函数 distHaversine() 转换为未向量化的 distm()。
如果您查看源代码，您会看到：
function (x, y, fun = distHaversine)
{
[...]
for (i in 1:n) {
dm[i, ] = fun(x[i, ], y)
}
return(dm)
}
虽然 distHaversine() 将“整个对象”发送到 C，但 distm() 将数据“按行”发送到 distHaversine()，因此强制 distHaversine() 执行在 C 中执行代码时也是如此。因此，不应使用 distm()。在性能方面，我看到使用包装函数 distm() 的危害更大，因为我看到了好处。
2。解释“解决方案”中的代码：
a) 分组：
您要分析每个组的数据：省份。
分组可以通过以下方式完成：split(df1, df1$province).
b) 对“柱块”进行分组
您想找到 lat/lon 的所有唯一组合。第一个猜测可能是 expand.grid()，但这不适用于多列。幸运的是，Flick 先生处理了这件事 。
然后你有一个 data.frame() 所有可能的组合，只需要使用
mean(distHaversine(...)).
我最初的想法是查看 distHaversine 的源代码并将其复制到我将与 proxy 一起使用的函数中。
那会像这样工作（注意 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++ 并添加了多线程作为奖励：
EDIT：原来 s2d 助手远非最佳，
此版本现在使用给定的公式 here.
EDIT2：我刚刚发现 RcppThread，
可用于检测用户中断。
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel,RcppThread)]]
#include <cstddef> // size_t
#include <math.h> // sin, cos, sqrt, atan2, pow
#include <vector>
#include <RcppThread.h>
#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++) {
if (RcppThread::isInterrupted(ind % static_cast<int>(1e5) == 0)) return;
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;
RcppParallel::parallelFor(0, size, hc, grain);
RcppThread::checkUserInterrupt();
return avg;
}
此代码不会分配任何中间矩阵，
它会简单地计算每对下三角的距离，并在最后累加平均值。
Kahan 求和部分请参见 here。
如果您将该代码保存在 haversine.cpp 中，
那么您可以执行以下操作：
library(dplyr)
library(Rcpp)
library(RcppParallel)
library(RcppThread)
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
如果您有数百万行，您可能需要等待很长时间...
我还应该提到，在通过 RcppParallel 创建的线程中检测到用户中断是不可能的，
所以如果你开始计算你应该等到它完成，
或完全重新启动 R/RStudio。
见上面的 EDIT2。
关于复杂性
根据你的实际数据和你的电脑有多少核，
您很可能最终要等待几天才能完成计算。
这个问题具有二次复杂性
（每个省，可以这么说）。
这一行：
int size = lon.length() * (lon.length() - 1) / 2;
表示必须执行的（半正弦）距离计算量。
因此，如果行数增加 n 倍，
粗略地说，计算次数增加了n^2 / 2倍。
没有办法优化这个；
如果不先实际计算每个数字，就无法计算 N 个数字的平均值，
你将很难找到比多线程 C++ 代码更快的东西，
所以你要么等待它结束，
或者在问题上投入更多的核心，
无论是单台机器还是多台机器一起工作。
否则无法解决这个问题。
鉴于您的数据有数百万行，这听起来像是一个 "XY" 问题。 IE。您真正需要的答案不是您提出的问题的答案。
我打个比方：如果你想知道森林中树木的平均高度，你不需要测量每一棵树。您只需测量一个足够大的样本，以确保您的估计值有足够高的概率接近您需要的真实平均值。
使用每栋房子与其他每栋房子的距离进行强力计算不仅会占用过多的资源（即使使用优化的代码），而且还会提供比您可能需要的更多的小数位，或者是由数据准确性证明（GPS 坐标通常最多只能在几米内正确）。
因此，我建议在样本量上进行计算，该样本量仅与您的问题要求的准确性水平所需的一样大。例如，以下内容将在几秒钟内提供对 200 万行的估计，这相当于 4 位有效数字。您可以通过增加样本大小来提高准确性，但考虑到 GPS 坐标本身的不确定性，我怀疑这样做是否合理。
sample.size=1e6
lapply(split(df1[3:4], df1$province),
function(x) {
s1 = x[sample(nrow(x), sample.size, T), ]
s2 = x[sample(nrow(x), sample.size, T), ]
mean(distHaversine(s1, s2))
})
一些要测试的大数据：
N=1e6
df1 <- data.frame(
province = c(rep(1,N),rep(2,N)),
house = 1:(2*N),
lat = c(rnorm(N,-76), rnorm(N,-85)),
lon = c(rnorm(N,39), rnorm(N,-55,2)))
要了解此方法的准确性，我们可以使用 bootstrapping。对于以下演示，我仅使用 100,000 行数据，以便我们可以在短时间内执行 1000 bootstrap 次迭代：
N=1e5
df1 <- data.frame(lat = rnorm(N,-76,0.1), lon = rnorm(N,39,0.1))
dist.f = function(i) {
s1 = df1[sample(N, replace = T), ]
s2 = df1[sample(N, replace = T), ]
mean(distHaversine(s1, s2))
}
boot.dist = sapply(1:1000, dist.f)
mean(boot.dist)
# [1] 17580.63
sd(boot.dist)
# [1] 29.39302
hist(boot.dist, 20)
即对于这些测试数据，平均距离为 17,580 +/- 29 m。这是 0.1% 的变异系数，对于大多数用途来说可能足够准确。正如我所说，如果确实需要，您可以通过增加样本量来提高准确性。
您可以使用矢量化版本的半正弦距离，例如：
dist_haversine_for_dfs <- function (df_x, df_y, lat, r = 6378137)
{
if(!all(c("lat", "lon") %in% names(df_x))) {
stop("parameter df_x does not have column 'lat' and 'lon'")
}
if(!all(c("lat", "lon") %in% names(df_y))) {
stop("parameter df_x does not have column 'lat' and 'lon'")
}
toRad <- pi/180
df_x <- df_x * toRad
df_y <- df_y * toRad
dLat <- df_y[["lat"]] - df_x[["lat"]]
dLon <- df_y[["lon"]] - df_x[["lon"]]
a <- sin(dLat/2) * sin(dLat/2) + cos(df_x[["lat"]]) * cos(df_y[["lat"]]) *
sin(dLon/2) * sin(dLon/2)
a <- pmin(a, 1)
dist <- 2 * atan2(sqrt(a), sqrt(1 - a)) * r
return(dist)
}
然后使用 data.table 和包 arrangements（为了更快地生成组合），您可以执行以下操作：
library(data.table)
dt <- data.table(df1)
ids <- dt[, {
comb_mat <- arrangements::combinations(x = house, k = 2)
list(house_x = comb_mat[, 1],
house_y = comb_mat[, 2])}, by = province]
jdt <- cbind(ids,
dt[ids$house_x, .(lon_x=lon, lat_x=lat)],
dt[ids$house_y, .(lon_y=lon, lat_y=lat)])
jdt[, dist := dist_haversine_for_dfs(df_x = jdt[, .(lon = lon.x, lat = lat.x)],
df_y = jdt[, .(lon = lon.y, lat = lat.y)])]
jdt[, .(mean_dist = mean(dist)), by = province]
输出
province mean_dist
1: 1 15379.21
2: 2 793612.04
我在下面添加了一个使用 spatialrisk 包的解决方案。这个包中的关键函数是用 C++ (Rcpp) 编写的，因此速度非常快。
library(data.table)
library(tidyverse)
library(spatialrisk)
library(optiRum)
# Expand grid
grid <- function(x){
df <- x[, lat, lon]
optiRum::CJ.dt(df, df)
}
由于输出的每个元素都是一个数据框，purrr::map_dfr 用于将它们行绑定在一起：
data.table(df1) %>%
split(.$province) %>%
map_dfr(grid, .id = "province") %>%
mutate(distm = spatialrisk::haversine(lat, lon, i.lat, i.lon)) %>%
filter(distm > 0) %>%
group_by(province) %>%
summarize(distm_mean = mean(distm))
输出：
province distm_mean
<chr> <dbl>
1 1 15379.
2 2 793612.
我想计算每个省多个房屋之间的平均地理距离。
假设我有以下数据。
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
如何计算省内所有房屋之间的距离并收集每个省的平均距离？
原始数据集每个省都有数百万个观测值，因此性能也是一个问题。
我的 10 美分。您可以：
# subset the province
df1 <- df1[which(df1$province==1),]
# get all combinations
all <- combn(df1$house, 2, FUN = NULL, simplify = TRUE)
# run your function and get distances for all combinations
distances <- c()
for(col in 1:ncol(all)) {
a <- all[1, col]
b <- all[2, col]
dist <- distm(c(df1$lon[a], df1$lat[a]), c(df1$lon[b], df1$lat[b]), fun = distHaversine)
distances <- c(distances, dist)
}
# calculate mean:
mean(distances)
# [1] 15379.21
这为您提供了该省的平均值，您可以将其与其他方法的结果进行比较。例如评论中提到的sapply：
df1 <- df1[which(df1$province==1),]
mean(sapply(split(df1, df1$province), dist))
# [1] 1.349036
如您所见，它给出了不同的结果，因为 dist 函数可以计算不同类型的距离（如欧几里得）而不能计算半正弦或其他 "geodesic" 距离。包 geodist 似乎有选项可以让你更接近 sapply:
library(geodist)
library(magrittr)
# defining the data
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))
# defining the function
give_distance <- function(resultofsplit){
distances <- c()
for (i in 1:length(resultofsplit)){
sdf <- resultofsplit
sdf <- sdf[[i]]
sdf <- sdf[c("lon", "lat", "province", "house")]
sdf2 <- as.matrix(sdf)
sdf3 <- geodist(x=sdf2, measure="haversine")
sdf4 <- unique(as.vector(sdf3))
sdf4 <- sdf4[sdf4 != 0] # this is to remove the 0-distances
mean_dist <- mean(sdf4)
distances <- c(distances, mean_dist)
}
return(distances)
}
split(df1, df1$province) %>% give_distance()
#[1] 15379.21 793612.04
例如该函数将为您提供每个省的平均距离值。现在，我没有设法让 give_distance 与 sapply 一起工作，但这应该已经更有效率了。
参考此
toCheck <- sapply(split(df1, df1$province), function(x){
combn(rownames(x), 2, simplify = FALSE)})
names(toCheck) <- sapply(toCheck, paste, collapse = " - ")
sapply(toCheck, function(x){
distm(df1[x[1],c("lon","lat")], df1[x[2],c("lon","lat")],
fun = distHaversine)
})
# 1 - 2 1 - 3 2 - 3 4 - 5 4 - 6 5 - 6
# 11429.10 22415.04 12293.48 634549.20 1188925.65 557361.28
如果每个省份的记录数相同，则此方法有效。如果不是这种情况，那么为 toCheck 分配适当名称的第二部分以及我们在最后如何使用它应该随着 toCheck 列表结构的变化而改变。它不关心数据集的顺序。
对于您的实际数据集，toCheck 将成为嵌套列表，因此您需要像下面这样调整函数；对于此解决方案，我还没有为 toCheck 命名。 （df2可以在答案末尾找到）。
df2 <- df2[order(df2$province),] #sorting may even improve performance
names(toCheck) <- paste("province", unique(df2$province))
toCheck <- sapply(split(df2, df2$province), function(x){
combn(rownames(x), 2, simplify = FALSE)})
sapply(toCheck, function(x){ sapply(x, function(y){
distm(df2[y[1],c("lon","lat")], df2[y[2],c("lon","lat")], fun = distHaversine)
})})
# $`province 1`
# [1] 11429.10 22415.04 1001964.84 12293.48 1013117.36 1024209.46
#
# $`province 2`
# [1] 634549.2 1188925.7 557361.3
#
# $`province 3`
# [1] 590083.2
#
# $`province 4`
# [1] 557361.28 547589.19 11163.92
您可以进一步获得每个省份的mean()。此外，如果需要，重命名嵌套列表的元素应该不难，这样您就可以知道每个距离对应于哪些房屋。
df2 <- data.frame(province = c(1, 1, 1, 2, 2, 2, 1, 3, 3, 4,4,4),
house = c(1, 2, 3, 4, 5, 6, 7, 10, 9, 8, 11, 12),
lat = c(-76.6, -76.5, -76.4, -75.4, -80.9, -85.7, -85.6, -76.4, -75.4, -80.9, -85.7, -85.6),
lon = c(39.2, 39.1, 39.3, 60.8, 53.3, 40.2, 40.1, 39.3, 60.8, 53.3, 40.2, 40.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() 取自
解释：
1.性能
我会避免使用 distm()，因为它将 vectorised 函数 distHaversine() 转换为未向量化的 distm()。
如果您查看源代码，您会看到：
function (x, y, fun = distHaversine)
{
[...]
for (i in 1:n) {
dm[i, ] = fun(x[i, ], y)
}
return(dm)
}
虽然 distHaversine() 将“整个对象”发送到 C，但 distm() 将数据“按行”发送到 distHaversine()，因此强制 distHaversine() 执行在 C 中执行代码时也是如此。因此，不应使用 distm()。在性能方面，我看到使用包装函数 distm() 的危害更大，因为我看到了好处。
2。解释“解决方案”中的代码：
a) 分组：
您要分析每个组的数据：省份。
分组可以通过以下方式完成：split(df1, df1$province).
b) 对“柱块”进行分组
您想找到 lat/lon 的所有唯一组合。第一个猜测可能是 expand.grid()，但这不适用于多列。幸运的是，Flick 先生处理了这件事
然后你有一个 data.frame() 所有可能的组合，只需要使用
mean(distHaversine(...)).
我最初的想法是查看 distHaversine 的源代码并将其复制到我将与 proxy 一起使用的函数中。
那会像这样工作（注意 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++ 并添加了多线程作为奖励：
EDIT：原来 s2d 助手远非最佳，
此版本现在使用给定的公式 here.
EDIT2：我刚刚发现 RcppThread， 可用于检测用户中断。
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel,RcppThread)]]
#include <cstddef> // size_t
#include <math.h> // sin, cos, sqrt, atan2, pow
#include <vector>
#include <RcppThread.h>
#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++) {
if (RcppThread::isInterrupted(ind % static_cast<int>(1e5) == 0)) return;
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;
RcppParallel::parallelFor(0, size, hc, grain);
RcppThread::checkUserInterrupt();
return avg;
}
此代码不会分配任何中间矩阵， 它会简单地计算每对下三角的距离，并在最后累加平均值。 Kahan 求和部分请参见 here。
如果您将该代码保存在 haversine.cpp 中，
那么您可以执行以下操作：
library(dplyr)
library(Rcpp)
library(RcppParallel)
library(RcppThread)
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
如果您有数百万行，您可能需要等待很长时间...
我还应该提到，在通过
见上面的 EDIT2。RcppParallel 创建的线程中检测到用户中断是不可能的，
所以如果你开始计算你应该等到它完成，
或完全重新启动 R/RStudio。
关于复杂性
根据你的实际数据和你的电脑有多少核， 您很可能最终要等待几天才能完成计算。 这个问题具有二次复杂性 （每个省，可以这么说）。 这一行：
int size = lon.length() * (lon.length() - 1) / 2;
表示必须执行的（半正弦）距离计算量。
因此，如果行数增加 n 倍，
粗略地说，计算次数增加了n^2 / 2倍。
没有办法优化这个；
如果不先实际计算每个数字，就无法计算 N 个数字的平均值，
你将很难找到比多线程 C++ 代码更快的东西，
所以你要么等待它结束，
或者在问题上投入更多的核心，
无论是单台机器还是多台机器一起工作。
否则无法解决这个问题。
鉴于您的数据有数百万行，这听起来像是一个 "XY" 问题。 IE。您真正需要的答案不是您提出的问题的答案。
我打个比方：如果你想知道森林中树木的平均高度，你不需要测量每一棵树。您只需测量一个足够大的样本，以确保您的估计值有足够高的概率接近您需要的真实平均值。
使用每栋房子与其他每栋房子的距离进行强力计算不仅会占用过多的资源（即使使用优化的代码），而且还会提供比您可能需要的更多的小数位，或者是由数据准确性证明（GPS 坐标通常最多只能在几米内正确）。
因此，我建议在样本量上进行计算，该样本量仅与您的问题要求的准确性水平所需的一样大。例如，以下内容将在几秒钟内提供对 200 万行的估计，这相当于 4 位有效数字。您可以通过增加样本大小来提高准确性，但考虑到 GPS 坐标本身的不确定性，我怀疑这样做是否合理。
sample.size=1e6
lapply(split(df1[3:4], df1$province),
function(x) {
s1 = x[sample(nrow(x), sample.size, T), ]
s2 = x[sample(nrow(x), sample.size, T), ]
mean(distHaversine(s1, s2))
})
一些要测试的大数据：
N=1e6
df1 <- data.frame(
province = c(rep(1,N),rep(2,N)),
house = 1:(2*N),
lat = c(rnorm(N,-76), rnorm(N,-85)),
lon = c(rnorm(N,39), rnorm(N,-55,2)))
要了解此方法的准确性，我们可以使用 bootstrapping。对于以下演示，我仅使用 100,000 行数据，以便我们可以在短时间内执行 1000 bootstrap 次迭代：
N=1e5
df1 <- data.frame(lat = rnorm(N,-76,0.1), lon = rnorm(N,39,0.1))
dist.f = function(i) {
s1 = df1[sample(N, replace = T), ]
s2 = df1[sample(N, replace = T), ]
mean(distHaversine(s1, s2))
}
boot.dist = sapply(1:1000, dist.f)
mean(boot.dist)
# [1] 17580.63
sd(boot.dist)
# [1] 29.39302
hist(boot.dist, 20)
即对于这些测试数据，平均距离为 17,580 +/- 29 m。这是 0.1% 的变异系数，对于大多数用途来说可能足够准确。正如我所说，如果确实需要，您可以通过增加样本量来提高准确性。
您可以使用矢量化版本的半正弦距离，例如：
dist_haversine_for_dfs <- function (df_x, df_y, lat, r = 6378137)
{
if(!all(c("lat", "lon") %in% names(df_x))) {
stop("parameter df_x does not have column 'lat' and 'lon'")
}
if(!all(c("lat", "lon") %in% names(df_y))) {
stop("parameter df_x does not have column 'lat' and 'lon'")
}
toRad <- pi/180
df_x <- df_x * toRad
df_y <- df_y * toRad
dLat <- df_y[["lat"]] - df_x[["lat"]]
dLon <- df_y[["lon"]] - df_x[["lon"]]
a <- sin(dLat/2) * sin(dLat/2) + cos(df_x[["lat"]]) * cos(df_y[["lat"]]) *
sin(dLon/2) * sin(dLon/2)
a <- pmin(a, 1)
dist <- 2 * atan2(sqrt(a), sqrt(1 - a)) * r
return(dist)
}
然后使用 data.table 和包 arrangements（为了更快地生成组合），您可以执行以下操作：
library(data.table)
dt <- data.table(df1)
ids <- dt[, {
comb_mat <- arrangements::combinations(x = house, k = 2)
list(house_x = comb_mat[, 1],
house_y = comb_mat[, 2])}, by = province]
jdt <- cbind(ids,
dt[ids$house_x, .(lon_x=lon, lat_x=lat)],
dt[ids$house_y, .(lon_y=lon, lat_y=lat)])
jdt[, dist := dist_haversine_for_dfs(df_x = jdt[, .(lon = lon.x, lat = lat.x)],
df_y = jdt[, .(lon = lon.y, lat = lat.y)])]
jdt[, .(mean_dist = mean(dist)), by = province]
输出
province mean_dist
1: 1 15379.21
2: 2 793612.04
我在下面添加了一个使用 spatialrisk 包的解决方案。这个包中的关键函数是用 C++ (Rcpp) 编写的，因此速度非常快。
library(data.table)
library(tidyverse)
library(spatialrisk)
library(optiRum)
# Expand grid
grid <- function(x){
df <- x[, lat, lon]
optiRum::CJ.dt(df, df)
}
由于输出的每个元素都是一个数据框，purrr::map_dfr 用于将它们行绑定在一起：
data.table(df1) %>%
split(.$province) %>%
map_dfr(grid, .id = "province") %>%
mutate(distm = spatialrisk::haversine(lat, lon, i.lat, i.lon)) %>%
filter(distm > 0) %>%
group_by(province) %>%
summarize(distm_mean = mean(distm))
输出：
province distm_mean
<chr> <dbl>
1 1 15379.
2 2 793612.