Znajdź wszystkich najbliższych sąsiadów w określonej odległości

Question

Mam dużą listę współrzędnych x i y, przechowywanych w tablicy numpy.

Coordinates = [[ 60037633 289492298] 
[ 60782468 289401668] 
[ 60057234 289419794]] 
... 
... 

Co chcę jest znalezienie wszystkich najbliższych sąsiadów w określonej odległości (powiedzmy 3 metry) i zapisać wynik tak, że później można zrobić kilka dalszych analiz na wynik.

W przypadku większości pakietów stwierdziłem, że należy zdecydować, ile NN należy znaleźć, ale chcę po prostu wszystko w ustalonej odległości.

Jak mogę osiągnąć coś takiego i jaki jest najszybszy i najlepszy sposób osiągnięcia czegoś takiego w przypadku dużego zestawu danych (kilka milionów punktów)?

Answer 1

można użyć scipy.spatial.cKDTree:

import numpy as np 
import scipy.spatial as spatial 
points = np.array([(1, 2), (3, 4), (4, 5)]) 
point_tree = spatial.cKDTree(points) 
# This finds the index of all points within distance 1 of [1.5,2.5]. 
print(point_tree.query_ball_point([1.5, 2.5], 1)) 
# [0] 

# This gives the point in the KDTree which is within 1 unit of [1.5, 2.5] 
print(point_tree.data[point_tree.query_ball_point([1.5, 2.5], 1)]) 
# [[1 2]] 

# More than one point is within 3 units of [1.5, 1.6]. 
print(point_tree.data[point_tree.query_ball_point([1.5, 1.6], 3)]) 
# [[1 2] 
# [3 4]] 

---

Oto przykład pokazujący, jak można znaleźć wszystkie najbliższych sąsiadów do tablicy punktów, z jednej rozmowy do point_tree.query_ball_point:

import numpy as np 
import scipy.spatial as spatial 
import matplotlib.pyplot as plt 
np.random.seed(2015) 

centers = [(1, 2), (3, 4), (4, 5)] 
points = np.concatenate([pt+np.random.random((10, 2))*0.5 
         for pt in centers]) 
point_tree = spatial.cKDTree(points) 

cmap = plt.get_cmap('copper') 
colors = cmap(np.linspace(0, 1, len(centers))) 
for center, group, color in zip(centers, point_tree.query_ball_point(centers, 0.5), colors): 
    cluster = point_tree.data[group] 
    x, y = cluster[:, 0], cluster[:, 1] 
    plt.scatter(x, y, c=color, s=200) 

plt.show()