Наибольшее расстояние между набором точек долготы/широты

Я думаю, что следующее может быть полезным приближением, которое масштабируется линейно, а не квадратично с количеством точек, и его довольно легко реализовать:

1. вычислить центр масс M точек
2. найти точку P₀, максимально удаленную от M
3. найти точку P₁, максимально удаленную от P₀
4. аппроксимировать максимальный диаметр расстоянием между P₀ и P₁

Это можно обобщить, повторив шаг 3 N раз и взяв расстояние между P_N-1 и P_N

Шаг 1 можно эффективно выполнить, аппроксимируя M как среднее значение долготы и широты, что нормально, когда расстояния «небольшие» и полюса достаточно далеко. Другие шаги можно было бы выполнить, используя формулу точного расстояния, но они намного быстрее, если координаты точек можно аппроксимировать лежащими на плоскости. Как только «удаленная пара» (надеюсь, пара с максимальным расстоянием) найдена, расстояние до нее можно пересчитать по точной формуле.

Примером аппроксимации может быть следующее: если φ(M) и λ(M) — широта и долгота центра масс, вычисляемые как Σφ(P)/n и Σλ(P)/n,

• x(P) = (λ(P) - λ(M) + C) cos(φ(P))
• y(P) = φ(P) - φ(M) [это только для ясности, также может быть просто y(P) = φ(P)]

где C обычно равно 0, но может быть ± 360 °, если набор точек пересекает линию λ = ± 180 °. Чтобы найти максимальное расстояние, вам просто нужно найти

• max((x(P_N) - x(P_N-1))² + (y(P_N sub>) - y(P_N-1))²)

(вам не нужен квадратный корень, потому что он монотонный)

То же самое преобразование координат можно использовать для повторения шага 1 (в новой системе координат), чтобы иметь лучшую начальную точку. Подозреваю, что при выполнении некоторых условий вышеописанные шаги (без повторения шага 3) всегда приводят к "истинной дальней паре" (моя терминология). Если бы я только знал, какие условия...

ИЗМЕНИТЬ:

Я ненавижу опираться на чужие решения, но кому-то придется.

По-прежнему сохраняя вышеуказанные 4 шага с необязательным (но, вероятно, полезным, в зависимости от типичного распределения баллов) повторением шага 3 и следуя решение Spacedman, выполнение вычислений в 3D преодолевает ограничения близости и расстояния от полюсов:

• х (Р) = грех (ф (Р))
• y(P) = cos(φ(P)) sin(λ(P))
• z(P) = cos(φ(P)) cos(λ(P))

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

Центр масс определяется выражением x(M) = Σx(P)/n и т. д., а максимальное значение, которое нужно искать, равно

• max((x(P_N) - x(P_N-1))² + (y(P_N sub>) - y(P_N-1))² + (z(P_N) - z(P_{N- 1}))²)

Итак: вы сначала преобразуете сферические координаты в декартовы, затем начинаете с центра масс, чтобы найти, по крайней мере, за два шага (шаги 2 и 3) самую удаленную точку от предыдущей точки. Вы можете повторять шаг 3 до тех пор, пока расстояние увеличивается, возможно, с максимальным количеством повторений, но это не уведет вас от локального максимума. Старт из центра масс тоже не очень поможет, если точки разбросаны по всей Земле.

РЕДАКТИРОВАТЬ 2:

Я достаточно выучил R, чтобы написать ядро ​​алгоритма (хороший язык для анализа данных!)

Для плоского приближения без учета задачи вокруг линии λ=±180°:

# input: lng, lat (vectors)
rad = pi / 180;
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i = which.max((x - mean(x))^2 + (y       )^2)
j = which.max((x - x[i]   )^2 + (y - y[i])^2)
# output: i, j (indices)

На моем ПК поиск индексов i и j для 1000000 точек занимает меньше секунды.
Следующая 3D-версия немного медленнее, но работает для любого распределения точек (и не нуждается в поправках, когда λ =±180° линия пересекается):

# input: lng, lat
rad = pi / 180
x = sin(lat * rad)
f = cos(lat * rad)
y = sin(lng * rad) * f
z = cos(lng * rad) * f
i = which.max((x - mean(x))^2 + (y - mean(y))^2 + (z - mean(z))^2)
j = which.max((x - x[i]   )^2 + (y - y[i]   )^2 + (z - z[i]   )^2)
k = which.max((x - x[j]   )^2 + (y - y[j]   )^2 + (z - z[j]   )^2) # optional
# output: j, k (or i, j)

Вычисление k может быть опущено (т. е. результат может быть получен с помощью i и j) в зависимости от данных и требований. С другой стороны, мои эксперименты показали, что дальнейший расчет индекса бесполезен.

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

ИЗМЕНИТЬ 3:

К сожалению, относительная погрешность плоскостной аппроксимации в крайних случаях может достигать 1-1/√3 ≅ 42,3%, что может быть неприемлемым, хотя и очень редким. Алгоритм можно модифицировать, чтобы иметь верхнюю границу примерно 20%, которую я получил с помощью циркуля и линейки (аналитическое решение громоздко). Модифицированный алгоритм находит пару точек с локально максимальным расстоянием, затем повторяет те же шаги, но на этот раз начиная с середины первой пары, возможно, находя другую пару:

# input: lng, lat
rad = pi / 180
x = (lng - mean(lng)) * cos(lat * rad)
y = (lat - mean(lat))
i.n_1 = 1 # n_1: n-1
x.n_1 = mean(x)
y.n_1 = 0 # = mean(y)
s.n_1 = 0 # s: square of distance
repeat {
   s = (x - x.n_1)^2 + (y - y.n_1)^2
   i.n = which.max(s)
   x.n = x[i.n]
   y.n = y[i.n]
   s.n = s[i.n]
   if (s.n <= s.n_1) break
   i.n_1 = i.n
   x.n_1 = x.n
   y.n_1 = y.n
   s.n_1 = s.n
}
i.m_1 = 1
x.m_1 = (x.n + x.n_1) / 2
y.m_1 = (y.n + y.n_1) / 2
s.m_1 = 0
m_ok  = TRUE
repeat {
   s = (x - x.m_1)^2 + (y - y.m_1)^2
   i.m = which.max(s)
   if (i.m == i.n || i.m == i.n_1) { m_ok = FALSE; break }
   x.m = x[i.m]
   y.m = y[i.m]
   s.m = s[i.m]
   if (s.m <= s.m_1) break
   i.m_1 = i.m
   x.m_1 = x.m
   y.m_1 = y.m
   s.m_1 = s.m
}
if (m_ok && s.m > s.n) {
   i = i.m
   j = i.m_1
} else {
   i = i.n
   j = i.n_1
}
# output: i, j

Алгоритм 3D можно изменить аналогичным образом. Можно (как в 2D, так и в 3D случае) начать еще раз с середины второй пары точек (если она найдена). Верхняя граница в этом случае "оставлена ​​в качестве упражнения для читателя" :-).

Сравнение модифицированного алгоритма с (слишком) простым алгоритмом показало для нормального и квадратного равномерного распределения почти удвоение времени обработки и уменьшение средней ошибки с 0,6% до 0,03% (порядок величины). . Дальнейший перезапуск со средней точки приводит к немного лучшей средней ошибке, но почти равной максимальной ошибке.

ИЗМЕНИТЬ 4:

Мне еще нужно изучить эту статью, но она похоже, что 20%, которые я нашел с помощью компаса и линейки, на самом деле составляют 1-1/√(5-2√3) ≅ 19,3%

Теорема № 1: порядок любых двух расстояний по большому кругу на поверхности земли такой же, как порядок расстояний по прямой линии между точками, где вы прокладываете туннель через землю.

Следовательно, превратите вашу широту в x, y, z на основе сферической земли произвольного радиуса или эллипсоида с заданными параметрами формы. Это пара синусов/косинусов на точку (не на пару точек).

Теперь у вас есть стандартная трехмерная задача, которая не зависит от вычисления гаверсинусов расстояний. Расстояние между точками просто евклидово (Пифагор в 3d). Нужен квадратный корень и несколько квадратов, и вы можете опустить квадратный корень, если вас интересуют только сравнения.

В этом могут помочь причудливые пространственные древовидные структуры данных. Или такие алгоритмы, как http://www.tcs.fudan.edu.cn/rudolf/Courses/Algorithms/Alg_ss_07w/Webprojects/Qinbo_diameter/2d_alg.htm (нажмите «Далее» для 3D-методов). Или код C++ здесь: http://valis.cs.uiuc.edu/~sariel/papers/00/diameter/diam_prog.html

Как только вы нашли свою пару максимальных расстояний, вы можете использовать формулу Хаверсина, чтобы получить расстояние вдоль поверхности для этой пары.

Вот наивный пример, который плохо масштабируется (как вы говорите), как вы говорите, но может помочь в построении решения на R.

## lonlat points
n <- 100
d <- cbind(runif(n, -180, 180), runif(n, -90, 90))


library(sp)
## distances on WGS84 ellipsoid
x <- spDists(d, longlat = TRUE)

## row, then column index of furthest points
ind <- c(row(x)[which.max(x)], col(x)[which.max(x)])

## maps
library(maptools)
data(wrld_simpl)
plot(as(wrld_simpl, "SpatialLines"), col = "grey")

points(d, pch = 16, cex = 0.5)

## draw the points and a line between  on the page
points(d[ind, ], pch = 16)
lines(d[ind, ], lwd = 2)


## for extra credit, draw the great circle on which the furthest points lie
library(geosphere)


lines(greatCircle(d[ind[1], ], d[ind[2], ]), col = "firebrick")

Пакет geosphere предоставляет больше возможностей для расчета расстояния, если это необходимо. См. ?spDists в sp для деталей, используемых здесь.

Вы не говорите нам, будут ли эти точки расположены в достаточно малой части земного шара. Для действительно глобальных наборов точек моим первым предположением будет запуск наивного алгоритма O (n ^ 2), возможно, повышение производительности за счет некоторой пространственной индексации (R *-деревья, восьмеричные деревья и т. д.). Идея состоит в том, чтобы заранее сгенерировать список n*(n-1) треугольников в матрице расстояний и по частям передать его в библиотеку быстрого расстояний, чтобы свести к минимуму ввод-вывод и оттока процессов. Haversine в порядке, вы также можете сделать это с помощью метода Винсенти (наибольший вклад во время выполнения вносит квадратичная сложность, а не (фиксированное количество) итераций в формуле Винсенти). В качестве примечания, на самом деле вам не нужен R для этого.

РЕДАКТИРОВАТЬ № 2: The Barequet-Har-Peled алгоритм (как указал Spacedman в своем ответе) имеет O((n +1/(e^3))log(1/e)) сложности для e>0, и его стоит изучить.

Для квазиплоской задачи это известно как «диаметр выпуклой оболочки» и состоит из трех частей:

1. Вычисление выпуклой оболочки с помощью развертки Грэма, равной O(n*log(n)) - на самом деле, следует попытаться преобразовать точки в поперечную проекцию Меркатора (используя центр тяжести точек в наборе данных).
2. Поиск противоположных точек с помощью алгоритма Rotating Calipers — линейный O(n) .
3. Нахождение наибольшего расстояния среди всех пар антиподов — линейный поиск, O(n).

Ссылка с псевдокодом и обсуждением: http://fredfsh.com/2013/05/03/convex-hull-and-its-diameter/

См. также обсуждение связанного вопроса здесь: -existing-points">https://gis.stackexchange.com/questions/17358/how-can-i-find-the-the-farthest-point-from-a-set-of-existing-points

EDIT: решение Spacedman указало мне на алгоритм Malandain-Boissonnat ( документ в формате pdf здесь ). Однако это хуже или так же, как наивный алгоритм O (n ^ 2) грубой силы.