Найти все координаты внутри круга в географических данных в Python

Question

У меня есть миллионы географических точек. Для каждого из них я хочу найти все «соседние точки», то есть все остальные точки в некотором радиусе, скажем, несколько сотен метров.

Существует наивное O (N ^ 2) решение этой проблемы - просто вычислите расстояние всех пар точек. Однако, поскольку я имею дело с надлежащей метрикой расстояния (географического расстояния), должен быть более быстрый способ сделать это.

Я хотел бы сделать это в Python. Одно решение, которое приходит на ум, - это использовать некоторую базу данных (mySQL с расширениями ГИС, PostGIS) и надеяться, что такая база данных позаботится об эффективном выполнении операции, описанной выше, с использованием некоторого индекса. Я бы предпочел что-то попроще, но мне не нужно создавать и изучать такие технологии.

Пара баллов

• Я буду выполнять операцию «найти соседей» миллионы раз
• Данные останутся статичными
• Поскольку проблема в некотором смысле проста, я хотел бы видеть, что они представляют собой код на python, который решает ее.

С точки зрения кода Python, я хочу что-то вроде:

points = [(lat1, long1), (lat2, long2) ... ] # this list contains millions lat/long tuples
points_index = magical_indexer(points)
neighbors = []
for point in points:
    point_neighbors = points_index.get_points_within(point, 200) # get all points within 200 meters of point
    neighbors.append(point_neighbors) 

Answer 1

Сообщенный Eamon, я придумал простое решение с использованием btrees, реализованное в SciPy.

from scipy.spatial import cKDTree
from scipy import inf

max_distance = 0.0001 # Assuming lats and longs are in decimal degrees, this corresponds to 11.1 meters
points = [(lat1, long1), (lat2, long2) ... ]
tree = cKDTree(points)

point_neighbors_list = [] # Put the neighbors of each point here

for point in points:
    distances, indices = tree.query(point, len(points), p=2, distance_upper_bound=max_distance)
    point_neighbors = []
    for index, distance in zip(indices, distances):
        if distance == inf:
            break
        point_neighbors.append(points[index])
    point_neighbors_list.append(point_neighbors)

Answer 2

scipy

Перво-наперво: существуют уже существующие алгоритмы для таких вещей, как дерево kd .Scipy имеет реализацию Python cKDtree , которая может найти все точки в данном диапазоне.

Бинарный поиск

В зависимости от того, что вы делаете, реализация чего-то подобного можетбыть нетривиальнымКроме того, создание дерева является довольно сложным (возможно, с большими накладными расходами), и вы можете избежать простого взлома, который я использовал ранее:

1. Вычислить PCA набора данных,Вы хотите повернуть набор данных так, чтобы наиболее значимое направление было первым, а ортогональное (менее большое) второе направление - ну, во-вторых, вторым.Вы можете пропустить это и просто выбрать X или Y, но это вычислительно дешево и обычно легко реализуемо.Если вы просто выберите X или Y, выберите направление с большей дисперсией.
2. Сортируйте точки по главному направлению (назовите это направление X).
3. Чтобы найти ближайшего соседа данногоpoint, найдите индекс точки, ближайшей к X, с помощью бинарного поиска (если точка уже находится в вашей коллекции, возможно, вы уже знаете этот индекс и вам не нужен поиск).Итеративно смотрите на следующую и предыдущую точки, сохраняя наилучшее совпадение и его расстояние от точки поиска.Вы можете перестать смотреть, когда разница в X больше или равна расстоянию до наилучшего совпадения (на практике обычно очень мало точек).
4. Чтобы найти все точки в данном диапазоне, выполнитетак же, как в шаге 3, за исключением того, что не останавливайтесь, пока разница в X не превысит диапазон.

По сути, вы выполняете предварительную обработку O (N log (N)), и для каждой точки примерно o(sqrt (N)) - или больше , если распределение ваших очков плохое.Если точки распределены примерно равномерно, число точек ближе к X, чем к ближайшему соседу, будет порядка квадратного корня из N. Это менее эффективно, если в пределах вашего диапазона много точек, но никогда не намного хуже, чем грубая сила.

Одним из преимуществ этого метода является то, что он все исполняемый при очень небольшом распределении памяти и в основном может быть выполнен с очень хорошей локализацией памяти, что означает, что он работает довольно хорошо, несмотря на очевидные ограничения.

Триангуляция Делоне

Другая идея: Триангуляция Делоне может сработать.Для триангуляции Делоне указано, что ближайший сосед любой точки является смежным узлом.Интуиция заключается в том, что во время поиска вы можете поддерживать кучу (приоритетную очередь) на основе абсолютного расстояния от точки запроса.Выберите ближайшую точку, убедитесь, что она находится в диапазоне, и, если это так, добавьте всех ее соседей.Я подозреваю , что невозможно пропустить какие-либо моменты, подобные этой, но вам нужно более внимательно изучить это, чтобы быть уверенным ...