【563】DBSCAN通过地理距离实现

　　提前计算好地理距离矩阵，然后将函数名复制到DBSCAN的函数里面。

import pandas as pd 
import numpy as np 
import folium
from scipy.spatial import ConvexHull
from math import radians, cos, sin, asin, sqrt
from sklearn.cluster import DBSCAN

def geodistance(lon1, lat1, lon2, lat2): # 经度1，纬度1，经度2，纬度2 （十进制度数）
    """
    大圆距离，great_circle
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # 将十进制度数转化为弧度
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
 
    # haversine公式
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # 地球平均半径，单位为公里
    return c * r * 1000    

def get_distance_matrix_from_array(points_array):
    """
    构建距离矩阵，每个点之间的 great_circle 距离
    """
    
    num = len(points_array)
    distance_matrix = np.zeros((num, num))
    for i in range(num):
        for j in range(num):
            if i == j:
                continue
            lng1, lat1 = points_array[i]
            lng2, lat2 = points_array[j]
            dis = geodistance(lng1, lat1, lng2, lat2)
            distance_matrix[i][j] = dis
    return distance_matrix

def DBSCAN_pts(aoi_points, eps, minpts):
    """
    minpts: 用数量作为计算
    aoi_points: 字典
    {'望京小区': [(0,1), (0,2), (1,2), (1,3)], ...}
    """
    for aoi, pts in aoi_points.items():
        distance_matrix = get_distance_matrix_from_array(pts)
        y_pred = DBSCAN(eps=eps, 
                      min_samples=minpts, 
                      metric='precomputed'
                     ).fit_predict(distance_matrix)

        # y_pred 的输出结果肯定有 -1，极端情况只有 -1
        # -1 的结果就是噪声点
        # 对于只有 -1 的情况，认为是点比较分散，无法聚类，是要删除掉的

        # 如果都为 -1，直接舍弃
        if len(set(y_pred.tolist())) == 1 and (y_pred[0] == -1):
            continue

        # 去掉 -1 的点
        tmp_pts = np.array(pts)[y_pred != -1]
        y_pred = y_pred[y_pred != -1]

        # 聚类点太少的簇可以删掉，初定为 5
        for i in range(y_pred.max() + 1):
            if list(y_pred).count(i) <= 5: 
                tmp_pts = tmp_pts[y_pred != i]
                y_pred = y_pred[y_pred != i] 

        aoi_points_dbscan_1[aoi] = tmp_pts.tolist()    
    
    return aoi_points_dbscan_1