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

Question

我有一组lng/lat坐标.计算集合中任意两点之间最大距离的有效方法是什么("最大直径",如果你愿意的话)？

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

编辑:这些点位于一个足够小的区域,用于测量携带移动设备的人在一天内活动的区域.

Answer 1

定理#1:沿着地球表面的任意两个大圆距离的排序与作为穿过地球的点之间的直线距离的排序相同.

因此,根据任意半径的球形地球或给定形状参数的椭球,将您的纬度转换为x,y,z.这是每点的几个正弦/余弦(不是每对点).

现在你有一个标准的3-d问题,不依赖于计算Haversine距离.点之间的距离只是欧几里德(在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 ₀和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,并遵循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和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情况下)从第二对点(如果找到)的中点再次重新开始.在这种情况下,上限是"留给读者的练习":-).

修正算法与(太)简单算法的比较表明,对于正态分布和方形均匀分布,处理时间几乎加倍,平均误差从.6%减少到.03%(数量级) .从中点进一步重新启动会导致平均误差稍微好一点,但最大误差几乎相等.

编辑4:

我还要研究这篇文章,但看起来我发现罗盘和直边的20%实际上是1-1 /√(5-2√3)≅19.3%