notebook.community

I have a 2D cartesian grid, non-equidistant. The coordinates of each cell are stored in X[M,N] en Y[M,N]. These could be lat/lon (spheric).



In [101]:

    
x = np.linspace(0,10, num = 500)
y = np.linspace(30,50, num=1000)
X, Y = np.meshgrid(x,y)
m, n = X.shape

We have a point p = (x,y)



In [102]:

    
p = (5.3, 33.2)

How do we find the index tuple (i,j) So that (X[i,j], Y[i,j]) is closer to p than other points in the grid [X,Y]

Manual distance computation



In [103]:

    
%%timeit
# compute distances for all points
dist = np.sqrt((X-p[0])**2 + (Y-p[1])**2) 
# lookup the minimum in the distance vector (without ravel works also)
idx = np.argmin(dist.ravel())
# lookup the i,j indices
i, j = np.unravel_index(idx, (m, n))









    



100 loops, best of 3: 7.82 ms per loop



In [104]:

    
X[i,j], Y[i,j]









    Out[104]:





(5.2905811623246493, 33.203203203203202)

Rtree



In [105]:

    
import rtree



In [106]:

    
index = rtree.Rtree()



In [107]:

    
for idx, (i, j) in enumerate(np.ndindex(m, n)):
    index.add(idx, (X[i,j], Y[i,j]))



In [108]:

    
%%timeit
idx = index.nearest(p).next()
i, j = np.unravel_index(item.id, (m, n))









    



10000 loops, best of 3: 168 µs per loop

KD tree



In [109]:

    
import scipy.spatial
index = scipy.spatial.cKDTree(data=np.c_[X.ravel(),Y.ravel()])



In [110]:

    
%%timeit
d, idx = index.query(p)
i, j = np.unravel_index(item.id, (m, n))









    



10000 loops, best of 3: 47.4 µs per loop



In [110]: