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

Posted 2023-03-05

【中文标题】一组经度/纬度点之间的最大距离

【英文标题】：Greatest distance between set of longitude/latitude points

【发布时间】：2013-05-31 20:16:59

【问题描述】：

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

一种幼稚的方法是使用Haversine formula计算每2个点之间的距离并得到最大值，但这显然不能很好地缩放。

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

【问题讨论】：

如果距离“小”（例如，几十英里/公里），更简单的公式将为解决方案提供更好的常数因子

你能举个例子吗？

几乎等同于***.com/questions/7129482/… 最近与最远应该是您的问题和那个问题之间的微不足道的区别。

查看 spDists 的 sp 包和有关距离计算的其他选项的 geosphere 包。

@hatchet：最接近与最远的区别并非微不足道

【参考方案1】：

定理 #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 方法）。或者这里的 C++ 代码：http://valis.cs.uiuc.edu/~sariel/papers/00/diameter/diam_prog.html

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

【讨论】：

正确，并且始终适用，而我的解决方案不是。唯一的缺点是它是 O(n log n)（但只有近似值才能比这更好）

虽然...在理论中，您的定理 #1 仅适用于完美的球形地球，而不适用于通用椭球...

您会注意到没有定理 #1 的证明 :) 我可能应该称其为未经证实的假设...我仍在尝试为椭圆体寻找反例...啊，对于扁球体，极距与赤道上截然相反的点之间的距离...

漂亮。要证明您的定理，请参考有关弦长和角度的公式（等效于球体上的大圆长）。 The relationship (on a unit circle) is "chord length = 2*sin(angle)"，从 0 到 pi 单调递增，这证明了你关于两个量的顺序相同的观点。

值得注意的是，引用的 C++ 代码 (valis.cs.uiuc.edu/~sariel/papers/00/diameter/diam_prog.html) 是最坏情况的二次函数，但据称“对于大多数实际输入”是线性的。一个不错的发现。

【参考方案2】：

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

计算点的质心 M 找到距离M最大的点P₀ 找到到P₀的最大距离的点P₁ 用 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，但如果点集穿过 λ=±180° 线，则可以为 ± 360°。要找到最大距离，您只需找到

max((x(P_N) - x(P_N-1))² + (y(P _N) - y(P_N-1))²)

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

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

编辑：

我讨厌建立在别人的解决方案上，但总得有人这样做。

仍然保持上述 4 个步骤，可选（但可能有益，取决于点的典型分布）重复第 3 步， 并关注solution of 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)

在我的 PC 上，查找索引 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：

我还得研究this article，但看起来我用指南针和直尺找到的 20% 实际上是 1-1/√(5-2√3) ≅ 19.3%

【讨论】：

也许你可以提供一个小的实际例子来说明这在r 中是如何工作的？ （OP 试图实现这一目标的语言）。

@SimonO101：对不起，我不知道 r :-(

@SimonO101：现在我知道了一点 R :-)

太棒了！很高兴听到它以及很好的解决方案和算法。来自我的 +1。

谢谢。快速简单的近似，非常适合我的问题。可以使用geosphere::distHaversine(c(lat[i], lng[i]), c(lat[j], lng[j]))

计算最终距离

【参考方案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 包在需要时提供了更多用于距离计算的选项。有关此处使用的详细信息，请参阅 sp 中的 ?spDists。

【讨论】：

+1 用于演示 sp 和 geosphere 机器。我觉得对于大量点，最快的搜索可能是：（1）将地球表面划分为网格； (2) 计算所有占用网格单元之间的最小和最大距离（使用它们的顶点）； (3) 只保留一组单元格中的点，这些点在整体上比任何其他单元格都更远；然后 (4) 对它们进行细分，重复步骤 2、3 和 4，直到点数被充分筛选。需要大量记账，但在大多数情况下应该运行得很快。

我在想类似的事情，你可以很容易地用光栅粗略一下，但今天不适合我。这是一个很好的问题，希望我有机会探索其中的一些想法（和沃尔特的）。我在 20000 点上尝试了这个，它可以通过，但它非常浪费，而且 50000 对于 16Gb RAM 来说太多了。 :)

【参考方案4】：

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

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

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

使用 Graham's scan 计算凸包，即 O(n*log(n)) - 事实上，应该尝试将点转换为横向墨卡托投影（使用数据集中点的质心）。 通过Rotating Calipers 算法查找对映点 - 线性 O(n)。 在所有对映对中查找最大距离 - 线性搜索，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 here 中的论文）。但是，这与蛮力朴素 O(n^2) 算法更差或相同。