Матрица расстояний с геосферой: избежать повторного исчисления

Question

Я хочу вычислить расстояние между всеми точками в очень большой матрице, используя distm из geosphere.

См. Минимальный пример:

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))

Можете ли вы помочь мне найтиболее эффективный способ вычислить все эти расстояния, избегая делать дважды каждый?

Answer 1

Если вы хотите вычислить все попарные расстояния для точек x, лучше использовать distm(x), а не distm(x,x).Функция distm возвращает одинаковую симметричную матрицу в обоих случаях, но когда вы передаете ей один аргумент, она знает, что матрица симметрична, поэтому она не будет выполнять ненужные вычисления.

Вы можете рассчитать время.

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 

Вы также можете посмотреть код R для geosphere::distm, чтобы проверить, что он обрабатывает два случая по-разному.

В стороне: быстрый поиск в Google находит parallelDist: вычисление матрицы параллельного расстоянияна КРАН.Геодезическое расстояние является опцией.

Answer 2

Использование combn() из базы R может быть немного проще и, вероятно, быстрее, чем загрузка дополнительных пакетов.Затем 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 

Взглядыхорошо.

Answer 3

Вы можете подготовить фрейм данных возможных комбинаций без повторений (с пакетами 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