Geosphere 距离矩阵:避免重复计算
Matrix of distances with Geosphere: avoid repeat calculus
我想使用 geosphere
中的 distm
来计算一个非常大的矩阵中所有点之间的距离。
看一个最小的例子:
library(geosphere)
library(data.table)
coords <- data.table(coordX=c(1,2,5,9), coordY=c(2,2,0,1))
distances <- distm(coords, coords, fun = distGeo)
问题是由于我正在计算的距离的性质,distm
返回一个对称矩阵,因此,我可以避免计算超过一半的距离:
structure(c(0, 111252.129800202, 497091.059564718, 897081.91986428,
111252.129800202, 0, 400487.621661164, 786770.053508848, 497091.059564718,
400487.621661164, 0, 458780.072878927, 897081.91986428, 786770.053508848,
458780.072878927, 0), .Dim = c(4L, 4L))
你能帮我找到一种更有效的方法来计算所有这些距离,避免每个距离做两次吗?
您可以准备一个没有重复的可能组合的数据框(使用 gtools
包)。然后计算这些对的距离。这是代码:
library(gtools)
library(geosphere)
library(data.table)
coords <- data.table(coordX = c(1, 2, 5, 9), coordY = c(2, 2, 0, 1))
pairs <- combinations(n = nrow(coords), r = 2, repeats.allowed = F, v = c(1:nrow(coords)))
distances <- apply(pairs, 1, function(x) {
distm(coords[x[1], ], coords[x[2], ], fun = distGeo)
})
# Construct distances matrix
dist_mat <- matrix(NA, nrow = nrow(coords), ncol = nrow(coords))
dist_mat[upper.tri(dist_mat)] <- distances
dist_mat[lower.tri(dist_mat)] <- distances
dist_mat[is.na(dist_mat)] <- 0
print(dist_mat)
结果:
[,1] [,2] [,3] [,4]
[1,] 0.0 111252.1 497091.1 400487.6
[2,] 111252.1 0.0 897081.9 786770.1
[3,] 497091.1 400487.6 0.0 458780.1
[4,] 897081.9 786770.1 458780.1 0.0
如果要计算点 x
的所有成对距离,最好使用 distm(x)
而不是 distm(x,x)
。 distm
函数 returns 在两种情况下都使用相同的对称矩阵,但是当您向它传递单个参数时,它知道该矩阵是对称的,因此它不会进行不必要的计算。
你可以计时。
library("geosphere")
n <- 500
xy <- matrix(runif(n*2, -90, 90), n, 2)
system.time( replicate(100, distm(xy, xy) ) )
# user system elapsed
# 61.44 0.23 62.79
system.time( replicate(100, distm(xy) ) )
# user system elapsed
# 36.27 0.39 38.05
您还可以查看 geosphere::distm
的 R 代码,检查它对这两种情况的处理方式是否不同。
旁白:快速 google 搜索发现 parallelDist
:CRAN 上的并行距离矩阵计算。测地距离是一个选项。
使用基础 R 中的 combn()
可能比加载额外的包更简单,而且可能更快。然后,distm()
使用 distGeo()
作为来源,所以使用后者应该会更快。
coords <- as.data.frame(coords) # this won't work with data.tables though
cbind(t(combn(1:4, 2)), unique(geosphere::distGeo(coords[combn(1:4, 2), ])))
# [,1] [,2] [,3]
# [1,] 1 2 111252.1
# [2,] 1 3 497091.1
# [3,] 1 4 897081.9
# [4,] 2 3 786770.1
# [5,] 2 4 400487.6
# [6,] 3 4 458780.1
我们可以用基准测试来检查它。
Unit: microseconds
expr min lq mean median uq max neval cld
distm 555.690 575.846 597.7672 582.352 596.1295 904.718 100 b
distGeo 426.335 434.372 450.0196 441.516 451.8490 609.524 100 a
看起来不错。
我想使用 geosphere
中的 distm
来计算一个非常大的矩阵中所有点之间的距离。
看一个最小的例子:
library(geosphere)
library(data.table)
coords <- data.table(coordX=c(1,2,5,9), coordY=c(2,2,0,1))
distances <- distm(coords, coords, fun = distGeo)
问题是由于我正在计算的距离的性质,distm
返回一个对称矩阵,因此,我可以避免计算超过一半的距离:
structure(c(0, 111252.129800202, 497091.059564718, 897081.91986428,
111252.129800202, 0, 400487.621661164, 786770.053508848, 497091.059564718,
400487.621661164, 0, 458780.072878927, 897081.91986428, 786770.053508848,
458780.072878927, 0), .Dim = c(4L, 4L))
你能帮我找到一种更有效的方法来计算所有这些距离,避免每个距离做两次吗?
您可以准备一个没有重复的可能组合的数据框(使用 gtools
包)。然后计算这些对的距离。这是代码:
library(gtools)
library(geosphere)
library(data.table)
coords <- data.table(coordX = c(1, 2, 5, 9), coordY = c(2, 2, 0, 1))
pairs <- combinations(n = nrow(coords), r = 2, repeats.allowed = F, v = c(1:nrow(coords)))
distances <- apply(pairs, 1, function(x) {
distm(coords[x[1], ], coords[x[2], ], fun = distGeo)
})
# Construct distances matrix
dist_mat <- matrix(NA, nrow = nrow(coords), ncol = nrow(coords))
dist_mat[upper.tri(dist_mat)] <- distances
dist_mat[lower.tri(dist_mat)] <- distances
dist_mat[is.na(dist_mat)] <- 0
print(dist_mat)
结果:
[,1] [,2] [,3] [,4]
[1,] 0.0 111252.1 497091.1 400487.6
[2,] 111252.1 0.0 897081.9 786770.1
[3,] 497091.1 400487.6 0.0 458780.1
[4,] 897081.9 786770.1 458780.1 0.0
如果要计算点 x
的所有成对距离,最好使用 distm(x)
而不是 distm(x,x)
。 distm
函数 returns 在两种情况下都使用相同的对称矩阵,但是当您向它传递单个参数时,它知道该矩阵是对称的,因此它不会进行不必要的计算。
你可以计时。
library("geosphere")
n <- 500
xy <- matrix(runif(n*2, -90, 90), n, 2)
system.time( replicate(100, distm(xy, xy) ) )
# user system elapsed
# 61.44 0.23 62.79
system.time( replicate(100, distm(xy) ) )
# user system elapsed
# 36.27 0.39 38.05
您还可以查看 geosphere::distm
的 R 代码,检查它对这两种情况的处理方式是否不同。
旁白:快速 google 搜索发现 parallelDist
:CRAN 上的并行距离矩阵计算。测地距离是一个选项。
使用基础 R 中的 combn()
可能比加载额外的包更简单,而且可能更快。然后,distm()
使用 distGeo()
作为来源,所以使用后者应该会更快。
coords <- as.data.frame(coords) # this won't work with data.tables though
cbind(t(combn(1:4, 2)), unique(geosphere::distGeo(coords[combn(1:4, 2), ])))
# [,1] [,2] [,3]
# [1,] 1 2 111252.1
# [2,] 1 3 497091.1
# [3,] 1 4 897081.9
# [4,] 2 3 786770.1
# [5,] 2 4 400487.6
# [6,] 3 4 458780.1
我们可以用基准测试来检查它。
Unit: microseconds
expr min lq mean median uq max neval cld
distm 555.690 575.846 597.7672 582.352 596.1295 904.718 100 b
distGeo 426.335 434.372 450.0196 441.516 451.8490 609.524 100 a
看起来不错。