在非矩形二维网格上高效地找到最近点的索引
Efficiently find indices of nearest points on non-rectangular 2D grid
我有一个不规则（非矩形）lon/lat 网格和 lon/lat 坐标中的一堆点，它们应该对应于网格上的点（尽管它们可能会稍微偏离数字原因）。现在我需要相应 lon/lat 点的索引。
我已经编写了一个函数来执行此操作，但它真的很慢。
def find_indices(lon,lat,x,y):
lonlat = np.dstack([lon,lat])
delta = np.abs(lonlat-[x,y])
ij_1d = np.linalg.norm(delta,axis=2).argmin()
i,j = np.unravel_index(ij_1d,lon.shape)
return i,j
ind = [find_indices(lon,lat,p*) for p in points]
我很确定 numpy/scipy 中有更好（更快）的解决方案。我已经用谷歌搜索了很多，但到目前为止我还没有找到答案。
有什么建议可以有效地找到相应（最近）点的索引吗？
PS: 这个问题来自另一个问题()。
编辑：解决方案
根据@Cong Ma的回答，我找到了以下解决方案：
def find_indices(points,lon,lat,tree=None):
if tree is None:
lon,lat = lon.T,lat.T
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
tree = sp.spatial.cKDTree(lonlat)
dist,idx = tree.query(points,k=1)
ind = np.column_stack(np.unravel_index(idx,lon.shape))
return [(i,j) for i,j in ind]
为了正确看待这个解决方案以及 Divakar 的回答中的解决方案，这里是我正在使用 find_indices 的函数的一些时间安排（以及它在速度方面的瓶颈所在）（请参阅link 以上）：
spatial_contour_frequency/pil0 : 331.9553
spatial_contour_frequency/pil1 : 104.5771
spatial_contour_frequency/pil2 : 2.3629
spatial_contour_frequency/pil3 : 0.3287
pil0 是我的初始方法，pil1 Divakar 的，pil2/pil3 是上面的最终解决方案，树是在运行中创建的pil2（即对于调用 find_indices 的循环的每次迭代）并且在 pil3 中仅一次（有关详细信息，请参阅其他线程）。尽管 Divakar 对我最初方法的改进使我的速度提高了 3 倍，但 cKDTree 将其提升到了一个全新的水平，速度又提高了 50 倍！将树的创建移出函数会使事情变得更快。
如果点足够本地化，你可以直接尝试scipy.spatial的cKDTree实现，正如我自己讨论的那样。 post 是关于插值的，但你可以忽略它并只使用查询部分。
tl;dr 版本：
阅读 scipy.sptial.cKDTree 的文档。通过将 (n, m) 形状的 numpy ndarray 对象传递给初始化程序来创建树，树将从 n m 维坐标创建。
tree = scipy.spatial.cKDTree(array_of_coordinates)
之后，使用 tree.query() 检索第 k 个最近的邻居（可能使用近似和并行化，请参阅文档），或使用 tree.query_ball_point() 查找给定范围内的所有邻居距离公差。
如果点没有很好地定位，并且出现球面曲率/非平凡拓扑，您可以尝试将流形分成多个部分，每个部分都小到可以被认为是局部的。
这是使用 scipy.spatial.distance.cdist -
的通用矢量化方法
import scipy
# Stack lon and lat arrays as columns to form a Nx2 array, where is N is grid**2
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
# Get the distances and get the argmin across the entire N length
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
# Get the indices corresponding to grid's shape as the final output
ind = np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
样本运行-
In [161]: lon
Out[161]:
array([[-11. , -7.82 , -4.52 , -1.18 , 2.19 ],
[-12. , -8.65 , -5.21 , -1.71 , 1.81 ],
[-13. , -9.53 , -5.94 , -2.29 , 1.41 ],
[-14.1 , -0.04 , -6.74 , -2.91 , 0.976]])
In [162]: lat
Out[162]:
array([[-11.2 , -7.82 , -4.51 , -1.18 , 2.19 ],
[-12. , -8.63 , -5.27 , -1.71 , 1.81 ],
[-13.2 , -9.52 , -5.96 , -2.29 , 1.41 ],
[-14.3 , -0.06 , -6.75 , -2.91 , 0.973]])
In [163]: lonlat = np.column_stack((lon.ravel(),lat.ravel()))
In [164]: idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
In [165]: np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
Out[165]: [[0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [3, 3]]
运行时测试 -
定义函数：
def find_indices(lon,lat,x,y):
lonlat = np.dstack([lon,lat])
delta = np.abs(lonlat-[x,y])
ij_1d = np.linalg.norm(delta,axis=2).argmin()
i,j = np.unravel_index(ij_1d,lon.shape)
return i,j
def loopy_app(lon,lat,pts):
return [find_indices(lon,lat,pts[i,0],pts[i,1]) for i in range(pts.shape[0])]
def vectorized_app(lon,lat,points):
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
return np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
时间安排：
In [179]: lon = np.random.rand(100,100)
In [180]: lat = np.random.rand(100,100)
In [181]: points = np.random.rand(50,2)
In [182]: %timeit loopy_app(lon,lat,points)
10 loops, best of 3: 47 ms per loop
In [183]: %timeit vectorized_app(lon,lat,points)
10 loops, best of 3: 16.6 ms per loop
为了获得更多性能，np.concatenate could be used in place of np.column_stack。
我有一个不规则（非矩形）lon/lat 网格和 lon/lat 坐标中的一堆点，它们应该对应于网格上的点（尽管它们可能会稍微偏离数字原因）。现在我需要相应 lon/lat 点的索引。
我已经编写了一个函数来执行此操作，但它真的很慢。
def find_indices(lon,lat,x,y):
lonlat = np.dstack([lon,lat])
delta = np.abs(lonlat-[x,y])
ij_1d = np.linalg.norm(delta,axis=2).argmin()
i,j = np.unravel_index(ij_1d,lon.shape)
return i,j
ind = [find_indices(lon,lat,p*) for p in points]
我很确定 numpy/scipy 中有更好（更快）的解决方案。我已经用谷歌搜索了很多，但到目前为止我还没有找到答案。
有什么建议可以有效地找到相应（最近）点的索引吗？
PS: 这个问题来自另一个问题(
编辑：解决方案
根据@Cong Ma的回答，我找到了以下解决方案：
def find_indices(points,lon,lat,tree=None):
if tree is None:
lon,lat = lon.T,lat.T
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
tree = sp.spatial.cKDTree(lonlat)
dist,idx = tree.query(points,k=1)
ind = np.column_stack(np.unravel_index(idx,lon.shape))
return [(i,j) for i,j in ind]
为了正确看待这个解决方案以及 Divakar 的回答中的解决方案，这里是我正在使用 find_indices 的函数的一些时间安排（以及它在速度方面的瓶颈所在）（请参阅link 以上）：
spatial_contour_frequency/pil0 : 331.9553
spatial_contour_frequency/pil1 : 104.5771
spatial_contour_frequency/pil2 : 2.3629
spatial_contour_frequency/pil3 : 0.3287
pil0 是我的初始方法，pil1 Divakar 的，pil2/pil3 是上面的最终解决方案，树是在运行中创建的pil2（即对于调用 find_indices 的循环的每次迭代）并且在 pil3 中仅一次（有关详细信息，请参阅其他线程）。尽管 Divakar 对我最初方法的改进使我的速度提高了 3 倍，但 cKDTree 将其提升到了一个全新的水平，速度又提高了 50 倍！将树的创建移出函数会使事情变得更快。
如果点足够本地化，你可以直接尝试scipy.spatial的cKDTree实现，正如我自己讨论的那样
tl;dr 版本：
阅读 scipy.sptial.cKDTree 的文档。通过将 (n, m) 形状的 numpy ndarray 对象传递给初始化程序来创建树，树将从 n m 维坐标创建。
tree = scipy.spatial.cKDTree(array_of_coordinates)
之后，使用 tree.query() 检索第 k 个最近的邻居（可能使用近似和并行化，请参阅文档），或使用 tree.query_ball_point() 查找给定范围内的所有邻居距离公差。
如果点没有很好地定位，并且出现球面曲率/非平凡拓扑，您可以尝试将流形分成多个部分，每个部分都小到可以被认为是局部的。
这是使用 scipy.spatial.distance.cdist -
import scipy
# Stack lon and lat arrays as columns to form a Nx2 array, where is N is grid**2
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
# Get the distances and get the argmin across the entire N length
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
# Get the indices corresponding to grid's shape as the final output
ind = np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
样本运行-
In [161]: lon
Out[161]:
array([[-11. , -7.82 , -4.52 , -1.18 , 2.19 ],
[-12. , -8.65 , -5.21 , -1.71 , 1.81 ],
[-13. , -9.53 , -5.94 , -2.29 , 1.41 ],
[-14.1 , -0.04 , -6.74 , -2.91 , 0.976]])
In [162]: lat
Out[162]:
array([[-11.2 , -7.82 , -4.51 , -1.18 , 2.19 ],
[-12. , -8.63 , -5.27 , -1.71 , 1.81 ],
[-13.2 , -9.52 , -5.96 , -2.29 , 1.41 ],
[-14.3 , -0.06 , -6.75 , -2.91 , 0.973]])
In [163]: lonlat = np.column_stack((lon.ravel(),lat.ravel()))
In [164]: idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
In [165]: np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
Out[165]: [[0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [0, 4], [3, 3]]
运行时测试 -
定义函数：
def find_indices(lon,lat,x,y):
lonlat = np.dstack([lon,lat])
delta = np.abs(lonlat-[x,y])
ij_1d = np.linalg.norm(delta,axis=2).argmin()
i,j = np.unravel_index(ij_1d,lon.shape)
return i,j
def loopy_app(lon,lat,pts):
return [find_indices(lon,lat,pts[i,0],pts[i,1]) for i in range(pts.shape[0])]
def vectorized_app(lon,lat,points):
lonlat = np.column_stack((lon.ravel(),lat.ravel()))
idx = scipy.spatial.distance.cdist(lonlat,points).argmin(0)
return np.column_stack((np.unravel_index(idx,lon.shape))).tolist()
时间安排：
In [179]: lon = np.random.rand(100,100)
In [180]: lat = np.random.rand(100,100)
In [181]: points = np.random.rand(50,2)
In [182]: %timeit loopy_app(lon,lat,points)
10 loops, best of 3: 47 ms per loop
In [183]: %timeit vectorized_app(lon,lat,points)
10 loops, best of 3: 16.6 ms per loop
为了获得更多性能，np.concatenate could be used in place of np.column_stack。