有效地找到非矩形2D网格上最近点的索引 [英] Efficiently find indices of nearest points on non-rectangular 2D grid

问题描述

我有一个不规则(非矩形)的lon/lat网格和lon/lat坐标中的一堆点，这些点应该对应于网格上的点(尽管由于数字原因它们可能会稍微偏离一点).现在，我需要相应的经度/纬度点的索引.

I have an irregular (non-rectangular) lon/lat grid and a bunch of points in lon/lat coordinates, which should correspond to points on the grid (though they might be slightly off for numerical reasons). Now I need the indices of the corresponding lon/lat points.

我编写了一个函数来执行此操作，但这确实很慢.

I've written a function which does this, but it is REALLY slow.

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中有更好(更快)的解决方案.我已经在Google上搜索了很多，但是到目前为止，答案还不知道我如何.

I'm pretty sure there's a better (and faster) solution in numpy/scipy. I've already googled quite a lot, but the answer has so far eluded me.

有人建议如何有效地找到相应(最近)点的索引吗?

PS:这个问题来自另一个问题(点击).

基于@Cong Ma的回答，我找到了以下解决方案:

Based on @Cong Ma's answer, I've found the following solution:

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的函数的一些时间安排(以及速度方面的瓶颈)(请参见上面的链接):

To put this solution and also the one from Divakar's answer into perspective, here are some timings of the function in which I'm using find_indices (and where it's the bottleneck in terms of speed) (see link above):

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倍，达到了一个全新的水平！将树的创建移出该函数可以使事情变得更快.

pil0 is my initial approach, pil1 Divakar's, and pil2/pil3 the final solution above, where the tree is created on-the-fly in pil2 (i.e. for every iteration of the loop in which find_indices is called) and only once in pil3 (see other thread for details). Even though Divakar's refinement of my initial approach gives me a 3x speed-up, cKDTree takes this to a whole new level with another 50x speedup! And moving the creation of the tree out of the function makes things even faster.

推荐答案

如果这些点已足够本地化，您可以直接尝试scipy.spatial的cKDTree实现，如我自己所讨论的

If the points are sufficiently localized, you may try directly scipy.spatial's cKDTree implementation, as discussed by myself in another post. That post was about interpolation but you can ignore that and just use the query part.

tl; dr版本:

阅读 scipy.sptial.cKDTree .通过将(n, m)形numpy ndarray对象传递给初始化程序来创建树，然后将根据n m维坐标创建树.

Read up the documentation of scipy.sptial.cKDTree. Create the tree by passing an (n, m)-shaped numpy ndarray object to the initializer, and the tree will be created from the n m-dimensional coordinates.

tree = scipy.spatial.cKDTree(array_of_coordinates)

然后，使用tree.query()检索第k个最近的邻居(可能具有近似和并行化功能，请参阅文档)，或使用tree.query_ball_point()查找在给定距离公差内的所有邻居.

After that, use tree.query() to retrieve the k-th nearest neighbor (possibly with approximation and parallelization, see docs), or use tree.query_ball_point() to find all neighbors within given distance tolerance.

如果这些点没有很好地定位，并且出现了球面曲率/非平凡的拓扑结构，则可以尝试将流形分成多个部分，每个部分足够小以至于被认为是局部的.

If the points are not well localized, and the spherical curvature / non-trivial topology kicks in, you can try breaking the manifold into multiple parts, each small enough to be considered local.

这篇关于有效地找到非矩形2D网格上最近点的索引的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！