r - 一组经度/纬度点之间的最大距离

Question

我有一组 lng/lat 坐标。计算集合中任意两点之间的最大距离（如果愿意的话，是“最大直径”）的有效方法是什么？

一种天真的方法是使用Haversine 公式计算每 2 个点之间的距离并获得最大值，但这显然不能很好地扩展。

编辑：这些点位于足够小的区域，测量携带移动设备的人在一天内活动的区域。

Answer 1

定理 #1：沿地球表面的任意两个大圆距离的排序与隧道穿过地球的点之间的直线距离的排序相同。

因此，根据任意半径的球形地球或给定形状参数的椭圆体，将您的 lat-long 转换为 x,y,z。这是每点的几个正弦/余弦（不是每对点）。

现在您有一个不依赖计算Haversine 距离的标准3-d 问题。点之间的距离只是欧几里得（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

找到最大距离对后，您可以使用 Haversine 公式获取该对沿表面的距离。

Answer 2

我认为以下可能是一个有用的近似值，它随点数线性而不是二次缩放，并且很容易实现：

1. 计算点的质心 M
2. 找到与M 有最大距离的点 P ₀
3. 找到到 P ₀距离最大的点 P ₁
4. _{以 P 0}和 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，但如果点集穿过 λ=±180° 线，则可以为 ± 360°。要找到最大距离，您只需找到

• 最大值((x(P _N ) - x(P _N-1 )) ² + (y(P _N ) - y(P _N-1 )) ² )

（你不需要平方根，因为它是单调的）

相同的坐标变换可用于重复步骤 1（在新坐标系中）以获得更好的起点。我怀疑如果满足某些条件，上述步骤（不重复步骤 3）总是会导致“真正的远距离对”（我的术语）。如果我只知道哪些条件...

编辑：

我讨厌建立在其他人的解决方案上，但有人必须这样做。

仍然保持上述 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 等给出，需要寻找的最大值是

• 最大值((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在我的 PC 上，查找索引和j1000000 点只需不到一秒钟的时间。
下面的 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%

Answer 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详细信息。

Answer 4

你没有告诉我们这些点是否会位于地球上足够小的一部分。对于真正的全局点集，我的第一个猜测是运行一个简单的 O(n^2) 算法，可能会通过一些空间索引（R*-trees、octal-trees 等）来提高性能。这个想法是在距离矩阵中预先生成一个 n*(n-1) 三角形列表，并将其以块的形式提供给快速距离库，以最大限度地减少 I/O 和进程流失。Haversine 很好，您也可以使用 Vincenty 的方法（运行时间的最大贡献者是二次复杂度，而不是 Vincenty 公式中的（固定数量的）迭代）。作为旁注，事实上，这些东西不需要 R。

编辑#2：Barequet-Har-Peled算法（正如 Spacedman 在他的回复中指出的那样）对于 e>0 具有 O((n+1/(e^3))log(1/e)) 复杂度，并且是值得探索。

对于准平面问题，这被称为“凸包直径”，它包含三个部分：

1. 使用Graham 的 O(n*log(n)) 扫描计算凸包- 事实上，应该尝试将点转换为横向墨卡托投影（使用数据集中点的质心）。
2. 通过旋转卡尺算法找到对映点- 线性 O(n)。
3. 在所有对映点对中找到最大距离 - 线性搜索，O(n)。

带有伪代码和讨论的链接：http: //fredfsh.com/2013/05/03/convex-hull-and-its-diameter/

另请参阅此处有关相关问题的讨论：https ://gis.stackexchange.com/questions/17358/how-can-i-find-the-farthest-point-from-a-set-of-existing-points

编辑：Spacedman 的解决方案向我指出了Malandain-Boissonnat算法（请参阅此处pdf 中的论文）。但是，这与蛮力朴素 O(n^2) 算法更差或相同。