现在有一个存有70个地址和城市名的文本，而没有这些地点的距离信息，而我们想要对这些地点进行聚类，找到每个簇的质心地点，从而可以安排合理的行程，即不同簇中的地点之间选择交通工具抵达，而位于同一个簇内的地点之间可以采取步行的方法抵达。使用Kmeans算法可以为我们找到一种更加经济而且高效的出行方式。

K-均值聚类算法实例——分类地图上的地点-最码农

数据说明：东经取经度的正值(Longitude)，西经取经度负值(-Longitude)，北纬取90-纬度值(90-Latitude),南纬取90+纬度值(90+Latitude)，则经过上述处理过后的两点被计为(MLonA, MLatA)和(MLonB, MLatB)。地点对应的经纬度经过以上处理过后，保存在places.txt文件和places.csv中的最后两列。

K-均值聚类算法实例——分类地图上的地点-最码农

使用球面余弦定理来计算两个经纬度之间的实际距离

地球平均半径6371千米。如果我们假设地球半径为R。设A点的经纬度为(LonA, LatA)，B点的经纬度为(LonB, LatB)，两点距离的如下公式：

  • a = sin(LatA*Pi/180)*sin(LatB*Pi/180)
  • b = cos(LatA*Pi/180)*cos(LatB*Pi/180)*cos((MLonA-MLonB)*Pi/180)
  • Distance = R* arccos(a + b)*Pi/180

# 计算根据经纬度计算两点之间的球面距离
def distSLC(vacA, vacB):
    a = sin(vacA[0, 1]*np.pi/180) * sin(vacB[0, 1]*np.pi/180)
    b = cos(vacA[0, 1]*np.pi/180) * cos(vacB[0, 1]*np.pi/180)*cos(np.pi*(vacB[0, 0]-vacA[0, 0])/180)
    return acos(a+b)*6371.0

K值的确定

本例使用拐点法，就是在不同的k值下计算簇内离差平方和，然后通过可视化的方法找到”拐点”所对应的k值，J为Kmeans算法的目标函数，随着簇数量的增加，簇中的样本量会越来越少，进而导致目标函数J的值也会越来越小，通过可视化方法，重点关注的是斜率的变化，当斜率由大突然变小时，并且之后的斜率变化缓慢，则认为突然变化的点就是寻找的目标点，因为继续随着簇数K的增加，聚类效果不再有大的变化。

完整代码及具体实现


import numpy as np
import pandas as pd
from math import cos, sin, acos
import matplotlib.pyplot as plt


# 加载数据
def loaddATEsET(filename):
    data = pd.read_csv(filename, header=None)
    dataArr = data.iloc[:, -2:].values
    return np.mat(dataArr)


# 计算根据经纬度计算两点之间的球面距离
def distSLC(vacA, vacB):
    a = sin(vacA[0, 1]*np.pi/180) * sin(vacB[0, 1]*np.pi/180)
    b = cos(vacA[0, 1]*np.pi/180) * cos(vacB[0, 1]*np.pi/180)*cos(np.pi*(vacB[0, 0]-vacA[0, 0])/180)
    return acos(a+b)*6371.0


# 随机初始化K个质心(质心满足数据边界之内)
def randCent(dataSet, k):
    n = np.shape(dataSet)[1]
    centroids = np.mat(np.zeros((k, n)))
    for j in range(n):
        minJ = np.min(dataSet[:, j])
        rangeJ = float(max(dataSet[:, j])-minJ)
        centroids[:, j] = minJ+rangeJ*np.random.rand(k, 1)
    return centroids


# Kmeans聚类算法
def kMeans(dataSet, k, distMeas=distSLC, createCent=randCent):
    m = np.shape(dataSet)[0]
    clusterAssment = np.mat(np.zeros((m, 2)))
    centroids = createCent(dataSet, k)
    clusterChanged = True
    while clusterChanged:
        clusterChanged = False
        for i in range(m):
            minDist = np.inf
            minIndex = -1
            for j in range(k):
                distJI = distMeas(centroids[j, :], dataSet[i, :])
                if distJI < minDist:
                    minDist = distJI
                    minIndex = j
            if clusterAssment[i, 0] != minIndex:
                clusterChanged = True
            clusterAssment[i, :] = minIndex, minDist**2
        # print(centroids)
        for cent in range(k):
            ptsInClust = dataSet[np.nonzero(clusterAssment[:, 0].A == cent)[0]]
            if len(ptsInClust) != 0:
                centroids[cent, :] = np.mean(ptsInClust, axis=0)
    return centroids, clusterAssment


# 计算Kmeans聚类算法的SSE值
def calcSSE(clusterAssment):
    return sum(clusterAssment[:, -1])


# 对地图坐标进行聚类，并在地图上显示聚类结果
def plotKMeans(dataSet, clusterAssment, centroids, k):
    fig = plt.figure()
    rect = [0.1, 0.1, 0.8, 0.8]
    scatterMarkers = ['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<']
    axprops = dict(xticks=[], yticks=[])
    ax0 = fig.add_axes(rect, label='ax0', **axprops)
    imgP = plt.imread('Portland.png')
    ax0.imshow(imgP)
    ax1 = fig.add_axes(rect, label='ax1', frameon=False)
    for i in range(k):
        ptsInCurrCluster = dataSet[np.nonzero(clusterAssment[:, 0].A == i)[0], :]
        markerStyle = scatterMarkers[i % len(scatterMarkers)]
        ax1.scatter(ptsInCurrCluster[:, 0].flatten().A[0],
                    ptsInCurrCluster[:, 1].flatten().A[0],
                    marker=markerStyle, s=90)
    ax1.scatter(centroids[:, 0].flatten().A[0], centroids[:, 1].flatten().A[0],
                marker='+', s=300)
    plt.show()
    # plt.savefig('fig.png', bbox_inches='tight')


dataSet = loaddATEsET("places.csv")
for i in range(7):
    centroids, clusterAssment = kMeans(dataSet, i+2)
    sseValue = calcSSE(clusterAssment)
print("当k值为：%d时，对应的SSE值为：%f" % (i+2, sseValue))
centroids, clusterAssment = kMeans(dataSet, 3)
plotKMeans(dataSet, clusterAssment, centroids, 3)

运行结果与分析

当k值分别为k=2,3,4,5,6,7,8时，计算出对应的SSE值如下：

K-均值聚类算法实例——分类地图上的地点-最码农
K取值2345678
SSE值3131.801819.441649.261077.58949.05762.87665.26

折线图如下：

K-均值聚类算法实例——分类地图上的地点-最码农

结论：

由上图可见，本例的K值确定为3。

效果图如下：

K-均值聚类算法实例——分类地图上的地点-最码农