距离大于零的Python KD树最近邻答案

【问题标题】：Python KD Tree Nearest Neigbour where distance is greater than zero距离大于零的Python KD树最近邻
【发布时间】：2014-07-20 11:41:24
【问题描述】：

我正在尝试对 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

【问题讨论】：

标签： python scipy spatial kdtree

【解决方案1】：

对于数据集中的每个点 P[i]，您都在问“我的数据集中哪个点最接近 P[i]？”你会得到答案“它是P[i]”。

如果你问一个不同的问题，“哪两个点最接近P[i]？”，即tree.query(pts,k=2)（与您的代码的区别是s/k=1/k=2/）你会得到P[i] 和P[j]，第二近的点，这就是你想要的结果。

旁注：

我建议您在构建树之前投影您的数据，因为在您的纬度范围内，经度 1 度距离的含义会有很大的波动。

【讨论】：

【解决方案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)

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

是 O(n^2)，在 10000 个点大约需要 1 秒
Ot 消耗大量内存（上述情况为 800 MB）

后一个问题当然可以通过分段来避免，但第一个问题不包括大点集。

这也可以使用scipy.spatial.distance.pdist：

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

这有点快（至少一半），但输出矩阵是压缩形式，请参阅scipy.spatial.distance.squareform 的文档。

如果您需要计算实际距离，那么这是一个不错的选择，因为pdist 可以处理球体上的距离。

然后，再次，您可以通过将查询扩展到两个最近点来使用您的 KDtree 方法：

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

然后nearest[1][0] 有点本身（“我、我自己和我”），nearest[1][1] 是真正最近的邻居（或者inf，如果没有足够近的地方）。

最佳解决方案取决于您拥有的点数。此外，如果您的地图点在地球上彼此不靠近，您可能希望使用笛卡尔 2D 距离以外的其他东西。

关于使用纬度和经度来查找距离的注意事项：如果您只是试图假装它们是二维笛卡尔点，那您就错了。在 60°N 处，一纬度为 1111 公里，而经度一度为 555 公里。因此，至少您必须将经度除以 cos（纬度）。即使有了这个技巧，当经度从东到西变化时，您最终也会遇到麻烦。

可能最简单的解决这个问题的方法是将坐标点计算为笛卡尔 3D 点：

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

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

【讨论】：