Python KD树最近邻，距离大于零

Question

我正在尝试为经纬度数据实现最近邻搜索。这是我的数据文件 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

Answer 1

怎么样，来个简单的解决方案？如果你有很多点（比如说一万个或更多），这种方法就不太合适了，但对于少量点来说，这种暴力解决方案可能会有用：

 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)，对于一万个点来说，大约需要一秒钟
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°N时，一度纬度是1111公里，而一度经度是555公里。所以，至少你需要把经度除以cos(纬度)。即使这样，当经度从东到西变化时，你也会遇到麻烦。

解决这个问题最简单的方法可能是将坐标点转换为三维笛卡尔坐标点：

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

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

Answer 2

对于你数据集中的每个点 P[i]，你在问“哪个点离 P[i] 最近？”然后你得到的答案是“就是 P[i] 自己”。

如果你问一个不同的问题，“离 P[i] 最近的两个点是什么？”也就是 tree.query(pts,k=2)（和你的代码的不同之处在于 s/k=1/k=2/），你会得到 P[i] 以及另一个点 P[j]，这个点是第二近的，这就是你想要的结果。

附带说明：

• 我建议你在构建树之前先对数据进行投影，因为在你的纬度范围内，1度经度的距离会有很大的波动。