按组划分的地理距离 - 在每对行上应用一个函数
问题描述 投票：3回答：6
我想计算每个省的一些房屋之间的平均地理距离。
假设我有以下数据。
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
如何计算省内所有房屋之间的距离并收集每个省的平均距离？
原始数据集每个省有数百万个观测值，因此性能也是一个问题。
6个回答
5
投票
投票
我最初的想法是查看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 ++并添加了多线程作为奖励：
编辑：原来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 ++代码更快的东西，所以你要么必须等待它，要么抛出更多的核心在这个问题上，无论是使用一台机器还是使用多台机器一起工作。否则你无法解决这个问题。
5
投票
投票
鉴于您的数据有数百万行，这听起来像是一个“XY”问题。即你真正需要的答案不是你问的问题的答案。
让我举个类比：如果你想知道森林中树木的平均高度，你就不会测量每棵树。您只需测量足够大的样本，以确保您的估计值足够高，可以根据需要接近真实平均值。
使用从每个房屋到每个其他房屋的距离来执行强力计算不仅会占用过多的资源（即使使用优化的代码），而且它将提供比您可能需要的更多的小数位，或者通过数据准确性来证明（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次自举迭代：
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％的变异系数，对于大多数目的而言可能足够准确。正如我所说，如果你真的需要，你可以通过增加样本量来获得更高的准确性。
4
投票
投票
解：
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.表现
我会避免使用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()，但这不适用于mulitple列。幸运的是，Flick先生照顾了这个expand.grid function for data.frames in R。
然后你有一个所有可能组合的data.frame()，只需要使用mean(distHaversine(...))。
1
投票
投票
关于这个thread，你问题的矢量化解决方案如下所示;
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))
0
投票
投票
我的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函数可以计算出不同类型的距离（如欧几里德）并且不能做半正弦或其他“测地”距离。包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合作，但这应该已经更有效了。
0
投票
投票
您可以使用漫游距离的矢量化版本，例如：
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
最新问题
- 是否可以使用Winscp连接到正在运行的docker容器？
- 为什么我的 Entra MSAL 流程适用于登录，但注销时却收到 404？
- Jetpack DataStore 本机库已添加到捆绑包中：libdatastore_shared_counter.so - 那是什么？
- Nextjs 将 prop 从服务器组件传递到客户端组件时出错
- DBT/Snowflake - 识别每个模型更新和插入的行
- 在包含空项的列表上使用带有谓词的 Exsist<T>
- 通过块将向量转换为矩阵 -reshape
- Git 日志 *.java 文件和提交者
- 需要通过Powershell安装Microsoft Visual C++ 2015-2022
- 包装类在我的 Eclipse 中不起作用，我错过了什么？
- pygame 中基于测验的游戏
- 配置 Gunicorn：未指定应用程序模块
- Python Click 模块，如何接受用户名和密码作为参数
- 在javascript中使用window.onbeforeunload事件中的window.event.keyCode捕获f5按键事件始终为0而不是116
- Pandas：获取特定数据类型的 value_count
- 如何输出与 HTML、CSS 和 JS 文件位于同一目录中的 Python 文件的结果
- 使用 Python 3 查找可用于 Gtk+3 小部件的信号/事件
- 自动对焦在 Material UI v5 中带有按钮组件的打开表单对话框中不起作用
- 如何在haskell中为优化编译器执行常量折叠算法？
- 使用 javascript 将 aria 标签添加到图像
© www.soinside.com 2019 - 2024. All rights reserved.