在r中获取大型数据集地理距离的最快方法

Question

我有两个包含经度和纬度的数据集。大数据集有大约20M的观测值，小数据集有36K的观测值。我试图从小数据集中找出落在大数据集中每个点200米范围内的观测数据的数量。所以过程是

1. 从大数据集中选取一个地理编码
2. 计算到小数据集中每个点的距离
3. 计算200米内有多少个点落在

问题是在r中完成这个任务需要很长时间，我使用geosphere包中的distm函数，但它仍然很慢。很容易超过24小时。

下面是我使用的代码:

library(geosphere)
#dataset1 (large)
df1 <- data.frame(longitude = c(-77.14239, -77.10750, -77.14239, -77.01797, -77.17203, -77.47230, -77.26490, -77.02824, -76.96993, -77.03185),
latitude = c(38.80575, 38.87987, 38.80575, 38.90425, 38.77076, 38.98140, 38.92800, 38.90436, 38.84569, 38.92080))
#dataset2 (small)
df2 <- data.frame(longitude = c(-75.34186, -123.59649, -108.20089, -115.16004, -87.62970),
latitude = c(40.11899, 44.38151, 36.71881, 36.22207, 41.71438))
# get the number of obs that fall within 200m
for(i in 1:nrow(df1)){
df1$n_within200m[i] <- sum(as.numeric(distm(cbind(df1$longitude[i], df1$latitude[i]),
cbind(df2$longitude, df2$latitude)) < 200))}

是否有更快的方法来完成这项任务?

Answer 1

这段代码给出了一个近似的答案，主要是为了解释这个想法。

首先，使用纬度和经度并不是存储数据的最佳方式，因为尺度与纬度不一致，靠近两极的情况会变得混乱。一个更好的方法是将数据存储为基于地球中心原点的直角坐标系。(这个坐标系通常被称为"以地球为中心的地球固定坐标系"或ECEF。)然后我们可以更简单地计算点之间的笛卡尔距离。我已经生成了我的测试数据直接在xyz坐标来节省时间-希望你不会有困难转换你的数据。

在大型数据集中搜索匹配项的一种快速方法是使用散列。我们将数据分成块，每边10公里，并将每个块的位置编码为唯一的整数。靠近的点通常会有相同的哈希值，对于任何给定的哈希值，在单个块中通常只有少量的对需要搜索。这个数字比搜索所有的数据要容易得多。

为了使解释简单，我只使用了一个散列。这意味着越过10公里街区边界的配对将被错过。为了解决这个问题，我们需要生成4个哈希，每个哈希在每个方向上偏移四分之一的块大小。H1 =散列(x, y, z), H2 =散列(2.5 x + y + 2.5, z + 2.5), H3 =散列(x + 5, y + 5, z + 5), H4 =散列(7.5 x + y + 7.5, z + 7.5)。如果两个点在任意哈希值上匹配，则需要检查对。生成每个散列的匹配列表，使用'rbind'组合列表，然后使用'unique'删除重复项。

Re <- 6371 # Earth's mean radius in km
# Define function to make some test data
# Uniformly distributed over surface of the Earth
# Include row number
# Lat/Long is inconvenient because distance is non-uniform around the poles
# So use Cartesian xyz instead
# Easy to convert lat/long to xyz
# Use (Re*cos(lat)*cos(lng),Re*cos(lat)*sin(lng),Re*sin(lat))
makeData <- function(n) {
pts <- data.frame(Row=seq(n),x=rnorm(n),y=rnorm(n),z=rnorm(n))
r <- sqrt(pts$x^2+pts$y^2+pts$z^2)
pts$x <- Re*pts$x/r
pts$y <- Re*pts$y/r
pts$z <- Re*pts$z/r
return(pts)
}
# Generate test data, 20 million points and 40 thousand points respectively
ptsA <- makeData(20000000)
ptsB <- makeData(40000)
# Generate an integer hash of the xyz position based on 10 km boxes
ptsA$hash <- floor(ptsA$x/10)+638*floor(ptsA$y/10)+638^2*floor(ptsA$z/10)
ptsB$hash <- floor(ptsB$x/10)+638*floor(ptsB$y/10)+638^2*floor(ptsB$z/10)
# Identify pairs of rows with a matching hash
matches <- merge(ptsA,ptsB,by="hash")
# Calculate distance between points for each pair identified
matches$dist <- with(matches,sqrt((x.x-x.y)^2+(y.x-y.y)^2+(z.x-z.y)^2))
# Filter out point pairs within 200 m of one another
matches <- matches[matches$dist < 0.2,]
# For each observed value of PtsA Row, count the number of matches
Counts <- aggregate(matches$Row.y,by=list(PtsARow=matches$Row.x),length)
# Initialise a column in PtsA with zeroes,
# then fill in the non-zero values found
ptsA$Count <- 0
ptsA[Counts$PtsARow,"Count"] <- Counts$x
# 'matches' contains a list of the matching points found
matches[order(matches$Row.x),][1:20,] # show first 20 matches
# Frequency table of counts
table(ptsA$Count)

我并不是说这是最快的方法，但它肯定会比24小时快得多。