Python KD树最近邻算法：距离大于零的最近邻。

Question

我正在尝试实现一个用于Lat和Lon数据的最近邻搜索。这是Data.txt文件。

61.3000183105 -21.2500038147 0
62.299987793 -23.750005722 1
66.3000488281 -28.7500038147 2
40.8000183105 -18.250005722 3
71.8000183105 -35.7500038147 3
39.3000183105 -19.7500019073 4
39.8000183105 -20.7500038147 5
41.3000183105 -20.7500038147 6

问题是，当我想对数据集上的每个经度和纬度执行最近邻操作时，它会搜索自身。例如，(-21.2500038147,61.3000183105) 的最近邻将是 (-21.2500038147,61.3000183105)，结果距离为0.0。我试图避免这种情况，但却没有成功。我尝试了 if not (array_equal)，但仍然无法解决问题。

以下是我的Python代码：

import numpy as np
from numpy import *
import decimal
from scipy import spatial
from scipy.spatial import KDTree
from math import radians,cos,sin,sqrt,exp


Lat =[]
Lon =[]
Day =[]

nja = []


Data = np.loadtxt('Data.txt',delimiter=" ")
for i in range(0,len(Data)):
    Lon.append(Data[i][:][0])
    Lat.append(Data[i][:][1])
    Day.append(Data[i][:][2])   

tree =spatial.KDTree(zip(Lon,Lat) )

print "Lon  :",len(Lon)
print "Tree :",len(tree.data)

for i in range(0,len(tree.data)):
    pts = np.array([tree.data[i][0],tree.data[i][1]])
    nja.append(pts)

for i in range(0, len(nja)):
    if not (np.array_equal(nja,tree.data)):
    nearest = tree.query(pts,k=1,distance_upper_bound =9)
    print nearest

- DBCooper

Answer 1

对于数据集中的每个点P [i]，您正在询问“在我的数据集中，哪个点最接近P [i]？”并且您会得到答案“它是P [i] ”。如果您提出不同的问题：“哪两个点最接近P [i]？”，即tree.query（pts，k = 2）（与您的代码不同之处在于s / k = 1 / k = 2 / ），您将获得P [i]和P [j]，即第二近的点，这就是您想要的结果。

附注：

我建议您在构建树之前对数据进行投影，因为在您的纬度范围内，经度距离的1度意味着变化很大。

- gboffi

Answer 2

如何使用低技术方案？如果您有大量的点（比如10000个或更多），这种方法就不再合理了，但对于较小数量的情况，这种暴力解决方案可能会有用：

 import numpy as np

 dist = (Lat[:,None]-Lat[None,:])**2 + (Lon[:,None]-Lon[None,:])**2

现在你有一个NxN的数组（N是点的数量），其中包含所有点对之间的距离（或距离的平方，更加精确）。然后，找到每个点的最短距离就是找到每行上的最小值。为了排除该点本身，您可以将对角线设置为NaN并使用nanargmax：

np.fill_diagonal(dist, np.nan)
closest = np.nanargmin(dist, axis=1)

这种方法非常简单，可以保证找到最接近的点，但有两个显著的缺点：

1. 它是O(n^2)的，在10000个点时需要大约一秒钟
2. 它消耗了大量的内存（对于上述情况，需要800 MB）

后一个问题当然可以通过分段处理来避免，但第一个问题排除了大型点集。

---

这也可以通过使用scipy.spatial.distance.pdist来实现：

dist=scipy.spatial.distance.pdist(np.column_stack((Lon, Lat)))

这会快一些（至少快了一半），但输出矩阵是在压缩形式下的，请查看 scipy.spatial.distance.squareform 的文档。

如果您需要计算真实距离，那么这是一个很好的选择，因为 pdist 可以处理球面上的距离。

---

那么，你可以通过将查询扩展到两个最近的点来再次使用你的KD树方法：

nearest = tree.query(pts, k=2, distance_upper_bound=9)

然后nearest[1][0]指的是自己这个点，nearest[1][1]则是真正的最近邻（如果没有足够近的话则为inf）。

最佳解决方案取决于您拥有的点的数量。此外，如果您的地图点彼此不接近，则可能希望使用其他东西而不是笛卡尔平面上的二维距离。

---

关于使用纬度和经度计算距离的注意事项：如果你只是把它们当作二维笛卡尔坐标点来处理，那么你会得到错误的结果。在北纬60°时，一度纬度相当于1111公里，而一度经度相当于555公里。因此，至少你需要将经度除以cos(纬度)。即使使用这个技巧，当经度从东向西变化时，你仍然会遇到麻烦。

也许最简单的解决方法是将坐标点计算为三维笛卡尔坐标点：

x = cos(lat) * cos(lon)
y = cos(lat) * sin(lon)
z = sin(lat)

如果您计算这些点之间的最短距离，您将得到正确的结果。（请注意，这些距离并不等同于地球表面上的真实最短距离。）

- DrV