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

Ответ 1

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

• вычислить центр масс M точек
• найдите точку P ₀, которая имеет максимальное расстояние до M
• найдите точку P ₁, которая имеет максимальное расстояние до P ₀
• приблизительный максимальный диаметр с расстоянием между P ₀ и P ₁

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

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

Примером аппроксимации может быть следующее: если φ (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) - y (P _N-1)) ²)

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

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

EDIT:

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

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

• x (P) = sin (φ (P))
• 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) - 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 точек.
Следующая трехмерная версия немного медленнее, но работает для любого распределения точек (и не необходимо изменить в случае пересечения линии λ = ± 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%

Ответ 2

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

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

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

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

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

Ответ 3

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

Ответ 4

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

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

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

• Вычислительная выпуклая оболочка с Graham scan, которая является O (n * log (n)) - на самом деле, нужно попробовать преобразовать точки в поперечная проекция Меркатора (с использованием центроида точек в наборе данных).
• Поиск антиподовых точек с помощью алгоритм вращающихся калибров - линейный O (n).
• Нахождение наибольшего расстояния между всеми антиподальными парами - линейный поиск, O (n).

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

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

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