查找两个大型数据集之间的最近坐标
Posted
技术标签:
【中文标题】查找两个大型数据集之间的最近坐标【英文标题】：Finding closest coordinates between two large data sets 【发布时间】：2019-09-09 03:32:14 【问题描述】：我的目标是根据两个数据集中的坐标确定数据集 2 中与数据集 1 中每个条目最近的条目。数据集 1 包含 180,000 行（只有 1,800 个唯一坐标），数据集 2 包含 4,500 行（完整的 4,500 个唯一坐标）。
我试图复制 *** 上类似问题的答案。例如：
R - Finding closest neighboring point and number of neighbors within a given radius, coordinates lat-long
Calculating the distance between points in different data frames
但是这些并不能以我想要的方式解决问题（它们要么加入数据帧，要么检查单个数据帧内的距离）。
Find the nearest X,Y coordinate using R 和 related posts 中的解决方案是我迄今为止找到的最接近的解决方案。
我对帖子的问题是它计算出单个数据帧内坐标之间的距离，我一直无法理解在 RANN::nn2 中更改哪些参数以跨两个数据帧进行。
建议的代码不起作用：
library(RANN)
dataset1[,4]<- nn2(data=dataset1, query=dataset2, k=2)
笔记/问题：
1) 应向查询提供哪个数据集以查找数据集 2 中与数据集 1 中给定值最接近的值？
2) 有什么方法可以避免数据集似乎需要相同宽度（列数）的问题？
3) 如何将输出（SRD_ID 和 distance）添加到数据集 1 的相关条目中？
4)RANN::nn2函数中eps参数有什么用？
目的是使用数据集 2 中最近的站 ID 以及数据集 1 中的条目与数据集 2 中的最近条目之间的距离填充数据集 1 中的 SRC_ID 和 distance 列。
下表展示了预期结果。 注意：SRC_ID 和 distance 值是我自己手动添加的示例值，几乎可以肯定是不正确的，并且可能不会被代码复制。
id HIGH_PRCN_LAT HIGH_PRCN_LON SRC_ID distance
1 3797987 52.88121 -2.873734 55 350
2 3798045 53.80945 -2.439163 76 2100
数据：
详情
platform x86_64-w64-mingw32
version.string R version 3.5.3 (2019-03-11)
数据集 1 输入（未缩小到唯一坐标）
structure(list(id = c(1L, 2L, 4L, 5L,
6L, 7L, 8L, 9, 10L, 3L),
HIGH_PRCN_LAT = c(52.881442267773, 57.8094538200198, 34.0233529,
63.8087900198, 53.6888144440184, 63.4462810678651, 21.6075544376207,
78.324442654172, 66.85532539759495, 51.623544596), HIGH_PRCN_LON = c(-2.87377812157822,
-2.23454414781635, -3.0984448341, -2.439163178635, -7.396111601421454,
-5.162345043546359, -8.63311254098095, 3.813289888829932,
-3.994325961186105, -8.9065532453272409), SRC_ID = c(NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA), distance = c(NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA)), row.names = c(NA, 10L), class = "data.frame")
数据集 2 输入
structure(list(SRC_ID = c(55L, 54L, 23L, 11L, 44L, 21L, 76L,
5688L, 440L, 61114L), HIGH_PRCN_LAT = c(68.46506, 50.34127, 61.16432,
42.57807, 52.29879, 68.52132, 87.83912, 55.67825, 29.74444, 34.33228
), HIGH_PRCN_LON = c(-5.0584, -5.95506, -5.75546, -5.47801, -3.42062,
-6.99441, -2.63457, -2.63057, -7.52216, -1.65532)), row.names = c(NA,
10L), class = "data.frame")
【问题讨论】：
如果他们是经纬度，您需要考虑使用geosphere::distHaversine 来获得正确的距离。但是按照答案中的建议找到最近的点。
【参考方案1】：
我写了一个引用这个thread 的答案。该函数被修改为负责报告距离并避免硬编码。请注意，它计算 欧几里得距离。
library(data.table)
#Euclidean distance
mydist <- function(a, b, df1, x, y)
dt <- data.table(sqrt((df1[[x]]-a)^2 + (df1[[y]]-b)^2))
return(data.table(Closest.V1 = which.min(dt$V1),
Distance = dt[which.min(dt$V1)]))
setDT(df1)[, j = mydist(HIGH_PRCN_LAT, HIGH_PRCN_LON, setDT(df2),
"HIGH_PRCN_LAT", "HIGH_PRCN_LON"),
by = list(id, HIGH_PRCN_LAT, HIGH_PRCN_LON)]
# id HIGH_PRCN_LAT HIGH_PRCN_LON Closest.V1 Distance.V1
# 1: 1 52.88144 -2.873778 5 0.7990743
# 2: 2 57.80945 -2.234544 8 2.1676868
# 3: 4 34.02335 -3.098445 10 1.4758202
# 4: 5 63.80879 -2.439163 3 4.2415854
# 5: 6 53.68881 -7.396112 2 3.6445416
# 6: 7 63.44628 -5.162345 3 2.3577811
# 7: 8 21.60755 -8.633113 9 8.2123762
# 8: 9 78.32444 3.813290 7 11.4936496
# 9: 10 66.85533 -3.994326 1 1.9296370
# 10: 3 51.62354 -8.906553 2 3.2180026
您可以使用RANN::nn2，但您需要确保使用正确的语法。以下作品！
as.data.frame(RANN::nn2(df2[,c(2,3)],df1[,c(2,3)],k=1))
# nn.idx nn.dists
# 1 5 0.7990743
# 2 8 2.1676868
# 3 10 1.4758202
# 4 3 4.2415854
# 5 2 3.6445416
# 6 3 2.3577811
# 7 9 8.2123762
# 8 7 11.4936496
# 9 1 1.9296370
# 10 2 3.2180026
【讨论】：
df1[ , c(4,5)] <- as.data.frame(RANN::nn2(df2[,c(2,3)],df1[,c(2,3)],k=1)) 将结果分配给所需的列。
感谢您的解决方案，尽管有数十万行，但使用 RANN::nn2 的第二个版本运行得非常快。您的评论中提供的代码段似乎将使用数据集 2 中最近条目的行号填充数据集 1 的第 4 列。您有什么建议可以改为输入数据集 2 中的 SRC_ID 值吗？跨度>
@Kickball：df1[ , c(4,5)] <- as.data.frame(RANN::nn2(df2[,c(2,3)],df1[,c(2,3)],k=1)); df1[,4] <- df2[df1[, 4], 1]
谢谢你，成功地将 SRC_ID 分配给数据集 1。我应该在前面提到它，但你知道是否有与 RANN:nn2 类似的函数，它使用理想的椭圆/文森特或球形/Haversine 距离?这对我的研究会更好。除了找到 geosphere::distVincentyEllipsoid 之外，我没有太多运气，但我不确定如何使代码适应不同的功能。如果您有任何想法，将不胜感激！
geosphere::distGeo 或相同的 raster::pointDistance 是比 Haversine 等更好（更精确）的选项。【参考方案2】：
数据
x = structure(list(id = c(1L, 2L, 4L, 5L,
6L, 7L, 8L, 9, 10L, 3L),
HIGH_PRCN_LAT = c(52.881442267773, 57.8094538200198, 34.0233529,
63.8087900198, 53.6888144440184, 63.4462810678651, 21.6075544376207,
78.324442654172, 66.85532539759495, 51.623544596), HIGH_PRCN_LON = c(-2.87377812157822,
-2.23454414781635, -3.0984448341, -2.439163178635, -7.396111601421454,
-5.162345043546359, -8.63311254098095, 3.813289888829932,
-3.994325961186105, -8.9065532453272409), SRC_ID = c(NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA), distance = c(NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA)), row.names = c(NA, 10L), class = "data.frame")
y = structure(list(SRC_ID = c(55L, 54L, 23L, 11L, 44L, 21L, 76L,
5688L, 440L, 61114L), HIGH_PRCN_LAT = c(68.46506, 50.34127, 61.16432,
42.57807, 52.29879, 68.52132, 87.83912, 55.67825, 29.74444, 34.33228
), HIGH_PRCN_LON = c(-5.0584, -5.95506, -5.75546, -5.47801, -3.42062,
-6.99441, -2.63457, -2.63057, -7.52216, -1.65532)), row.names = c(NA,
10L), class = "data.frame")
解决方案。请注意“3:2”以获取“经度/纬度”的顺序。
library(raster)
d <- pointDistance(x[,3:2], y[,3:2], lonlat=TRUE, allpairs=T)
i <- apply(d, 1, which.min)
x$SRC_ID = y$SRC_ID[i]
x$distance = d[cbind(1:nrow(d), i)]
x
# id HIGH_PRCN_LAT HIGH_PRCN_LON SRC_ID distance
#1 1 52.88144 -2.873778 44 74680.48
#2 2 57.80945 -2.234544 5688 238553.51
#3 4 34.02335 -3.098445 61114 137385.18
#4 5 63.80879 -2.439163 23 340642.70
#5 6 53.68881 -7.396112 44 308458.73
#6 7 63.44628 -5.162345 23 256176.88
#7 8 21.60755 -8.633113 440 908292.28
#8 9 78.32444 3.813290 76 1064419.47
#9 10 66.85533 -3.994326 55 185119.29
#10 3 51.62354 -8.906553 54 251580.45
图解
plot(x[,3:2], ylim=c(0,90), col="blue", pch=20)
points(y[,3:2], col="red", pch=20)
for (i in 1:nrow(x))
j <- y$SRC_ID==x$SRC_ID[i]
arrows(x[i,3], x[i,2], y[j,3], y[j,2],length=.1)
text(x[,3:2], labels=x$id, pos=1, cex=.75)
text(y[,3:2], labels=y$SRC_ID, pos=3, cex=.75)
【讨论】：
感谢您的详细回复，包括可视化！不幸的是，当我针对完整大小的数据集运行您的代码时，在i <- apply(d, 1, which.min) 阶段引发以下错误，Error: cannot allocate vector of size 5.7 Gb 由于 r 消耗了超过 11GB 的系统 RAM。我可以将我的唯一坐标对拆分为单独的数据框/集，然后将其用作数据集 1，但我会先尝试其他解决方案。以上是关于查找两个大型数据集之间的最近坐标的主要内容，如果未能解决你的问题，请参考以下文章