как мне сгруппировать список географических точек по расстоянию?

Question

У меня есть список точек P = [p1, ... pN], где pi = (широта I, долгота I).

Используя Python 3, я хотел бы найти наименьший набор кластеров (непересекающихся подмножеств P), чтобы каждый член кластера находился в пределах 20 км от любого другого члена кластера.

Расстояние между двумя точками рассчитывается по методу Винсенти .

Чтобы сделать это немного более конкретным, предположим, у меня есть набор точек, таких как

from numpy import *
points = array([[33.    , 41.    ],
       [33.9693, 41.3923],
       [33.6074, 41.277 ],
       [34.4823, 41.919 ],
       [34.3702, 41.1424],
       [34.3931, 41.078 ],
       [34.2377, 41.0576],
       [34.2395, 41.0211],
       [34.4443, 41.3499],
       [34.3812, 40.9793]])

Тогда я пытаюсь определить эту функцию:

from geopy.distance import vincenty
def clusters(points, distance):
    """Returns smallest list of clusters [C1,C2...Cn] such that
       for x,y in Ci, vincenty(x,y).km <= distance """
    return [points]  # Incorrect but gives the form of the output

ПРИМЕЧАНИЕ. Многие вопросы группируются по географическому местоположению и атрибуту . Мой вопрос только для местоположения . Это для широты / долготы, , а не Евклидово расстояние. Есть и другие вопросы, которые дают своего рода ответы, но не дают ответа на этот вопрос (многие без ответа):

• https://datascience.stackexchange.com/questions/761/clustering-geo-location-coordinates-lat-long-pairs
• https://gis.stackexchange.com/questions/300171/clustering-geo-points-and-export-borders-in-kml
• https://gis.stackexchange.com/questions/194873/clustering-geographical-data-based-on-point-location-and-associated-point-values
• https://gis.stackexchange.com/questions/256477/clustering-latitude-longitude-data-based-on-distance
• и более, ни один из которых не отвечает на этот вопрос.

Answer 1

почему бы не использовать библиотеку S2 для создания 20-километровых зон и посмотреть, какие точки находятся в каждой?

Answer 2

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

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

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

import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
import math

points = np.array([[33.    , 41.    ],
       [33.9693, 41.3923],
       [33.6074, 41.277 ],
       [34.4823, 41.919 ],
       [34.3702, 41.1424],
       [34.3931, 41.078 ],
       [34.2377, 41.0576],
       [34.2395, 41.0211],
       [34.4443, 41.3499],
       [34.3812, 40.9793]])


def distance(origin, destination): #found here https://gist.github.com/rochacbruno/2883505
    lat1, lon1 = origin[0],origin[1]
    lat2, lon2 = destination[0],destination[1]
    radius = 6371 # km
    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c

    return d

def create_clusters(number_of_clusters,points):
    kmeans = KMeans(n_clusters=number_of_clusters, random_state=0).fit(points)
    l_array = np.array([[label] for label in kmeans.labels_])
    clusters = np.append(points,l_array,axis=1)
    return clusters

def validate_solution(max_dist,clusters):
    _, __, n_clust = clusters.max(axis=0)
    n_clust = int(n_clust)
    for i in range(n_clust):
        two_d_cluster=clusters[clusters[:,2] == i][:,np.array([True, True, False])]
        if not validate_cluster(max_dist,two_d_cluster):
            return False
        else:
            continue
    return True

def validate_cluster(max_dist,cluster):
    distances = cdist(cluster,cluster, lambda ori,des: int(round(distance(ori,des))))
    print(distances)
    print(30*'-')
    for item in distances.flatten():
        if item > max_dist:
            return False
    return True

if __name__ == '__main__':
    for i in range(2,len(points)):
        print(i)
        print(validate_solution(20,create_clusters(i,points)))

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

Вы можете заменить лямбда-функцию в cdist на любую метрику расстояния, которую вы выбрали, я нашел большое расстояние в круге в упомянутом репо.

Answer 3

Вот решение, которое кажется правильным и будет вести себя в худшем случае O (N ^ 2) и лучше в зависимости от данных:

def my_cluster(S,distance):
    coords=set(S)
    C=[]
    while len(coords):
        locus=coords.pop()
        cluster = [x for x in coords if vincenty(locus,x).km <= distance]
        C.append(cluster+[locus])
        for x in cluster:
            coords.remove(x)
    return C

ПРИМЕЧАНИЕ : я не помечаю это как ответ, потому что одно из моих требований - это наименьший набор кластеров. Мой первый проход хорош, но я не доказал, что это самый маленький сет.

Результат (по большему набору точек) можно визуализировать следующим образом: