3. 對地圖上的點進行聚類

示例：對于地理數據應用二分K-均值算法

(1)收集數據：使用Yahoo! placeFinder API收集數據。
(2)準備數據：只保留經緯度信息。
(3)分析數據：使用matplotlib來構建一個二維數據圖，其中包含簇與地圖。
(4)訓練算法：訓練不適用無監督學習。
(5)測試算法：使用上篇中的bikmeans（）函教。
(6)使用算法：最后的輸出是包含簇及簇中心的地圖。

3.1 收集數據

# -*- coding: utf-8 -*- import urllib import json def geoGrab(stAddress, city): apiStem = 'http://where.yahooapis.com/geocode?' params = {}params['flags'] = 'J' # 返回JSONparams['appid'] = 'aaa0VN6k'params['location'] = '%s %s' % (stAddress, city)url_params = urllib.urlencode(params) # 將params字典轉換為可以通過URL進行傳遞的字符串格式yahooApi = apiStem + url_params print yahooApi # 輸出URLc=urllib.urlopen(yahooApi) #讀取返回值return json.loads(c.read()) # 返回一個字典from time import sleep def massPlaceFind(fileName):fw = open('places.txt', 'w')for line in open(fileName).readlines():line = line.strip()lineArr = line.split('\t') # 是以tab分隔的文本文件retDict = geoGrab(lineArr[1], lineArr[2]) # 讀取2列和第3列if retDict['ResultSet']['Error'] == 0: # 檢查輸出字典，判斷有沒有出錯lat = float(retDict['ResultSet']['Results'][0]['latitude']) # 讀取經緯度lng = float(retDict['ResultSet']['Results'][0]['longitude'])print "%s\t%f\t%f" % (lineArr[0], lat, lng)fw.write('%s\t%f\t%f\n' % (line, lat, lng)) # 添加到對應的行上else: print 'error fetching' # 有錯誤時不需要抽取經緯度sleep(1) # 避免頻繁調用API，過于頻繁的話請求會被封掉fw.close()geoGrab('1 VA Center', 'Augusta,ME')

Yahoo! placeFinder API：

為了使用該服務，需要注冊以獲得一個APIkey。具體地，需要在Yahoo！開發者網絡（Http://developer.yahoo.com)中進行注冊。創建一個桌面應用后會獲得一個appid。需要appid來使用geocoder。一個geocoder接受給定地址，然后返回該地址對應的經緯度。

JSON：

是一種用于序列化數組和字典的文件格式。JSON是javascript object
Notation的縮寫，接下來使用urllib的urlencode（）函數將創建的字典轉換為可以通過URL進行傳遞的字符串格式。最后，打開URL讀取返回值。由于返回值是json格式的，所以可以使用JSON的Python模塊來將其解碼為一個字典。一旦返回了解碼后的字典，也就意味著成功地對一個地址進行了地理編碼。

這里沒能實現，只是對代碼分析了下。。

3.2 對地理坐標進行聚類

# -*- coding: utf-8 -*- from numpy import * # K-均值聚類支持函數 def loadDataSet(fileName): dataMat = [] fr = open(fileName)for line in fr.readlines():curLine = line.strip().split('\t')fltLine = map(float,curLine) dataMat.append(fltLine)return dataMat# 計算兩個向量的歐式距離 def distEclud(vecA, vecB):return sqrt(sum(power(vecA - vecB, 2))) # 為給定數據集構建一個包含k個隨機質心的集合,是以每列的形式生成的 def randCent(dataSet, k):n = shape(dataSet)[1]centroids = mat(zeros((k,n))) for j in range(n): minJ = min(dataSet[:,j]) # 找到每一維的最小rangeJ = float(max(dataSet[:,j]) - minJ) # 每一維的最大和最小值之差centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1)) # 生成隨機值#print centroids[:,j]return centroids # 返回隨機質心,是和數據點相同的結構# k--均值聚類算法(計算質心--分配--重新計算) def kMeans(dataSet, k, distMeas=distEclud, createCent=randCent): # k是簇的數目m = shape(dataSet)[0] # 得到樣本的數目clusterAssment = mat(zeros((m,2))) # 創建矩陣來存儲每個點的簇分配結果# 第一列：記錄簇索引值，第二列：存儲誤差，歐式距離的平方centroids = createCent(dataSet, k) # 創建k個隨機質心clusterChanged = Truewhile clusterChanged: # 迭代使用while循環來實現clusterChanged = False for i in range(m): # 遍歷每個數據點，找到距離每個點最近的質心minDist = inf; minIndex = -1for j in range(k): # 尋找最近的質心distJI = distMeas(centroids[j,:],dataSet[i,:])if distJI < minDist:minDist = distJI; minIndex = jif clusterAssment[i,0] != minIndex: # 更新停止的條件clusterChanged = TrueclusterAssment[i,:] = minIndex,minDist**2 # minDist**2就去掉了根號 for cent in range(k): # 更新質心的位置ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]] centroids[cent,:] = mean(ptsInClust, axis=0) # 然后計算均值，axis=0:沿列方向 #print 'centroids:',centroidsreturn centroids, clusterAssment # 返回簇和每個簇的誤差值，誤差值是當前點到該簇的質心的距離# 二分k--均值聚類算法 def biKmeans(dataSet, k, distMeas=distEclud):m = shape(dataSet)[0] clusterAssment = mat(zeros((m,2))) # 存儲數據集中每個點的簇分配結果及平方誤差centroid0 = mean(dataSet, axis=0).tolist()[0] # 計算整個數據集的質心：1*2的向量centList =[centroid0] # []的意思是使用一個列表保存所有的質心,簇列表,[]的作用很大for j in range(m): # 遍歷所有的數據點，計算到初始質心的誤差值，存儲在第1列clusterAssment[j,1] = distMeas(mat(centroid0), dataSet[j,:])**2while (len(centList) < k): # 不斷對簇進行劃分，直到klowestSSE = inf # 初始化SSE為無窮大for i in range(len(centList)): # 遍歷每一個簇#print 'i:',i # 數組過濾得到所有的類別簇等于i的數據集ptsInCurrCluster = dataSet[nonzero(clusterAssment[:,0].A==i)[0],:]# 得到2個簇和每個簇的誤差，centroidMat：簇矩陣 splitClustAss：[索引值,誤差]centroidMat, splitClustAss = kMeans(ptsInCurrCluster, 2, distMeas) # centroidMat是矩陣sseSplit = sum(splitClustAss[:,1]) # 求二分k劃分后所有數據點的誤差和 # 數組過濾得到整個數據點集的簇中不等于i的點集#print nonzero(clusterAssment[:,0].A!=i)[0]sseNotSplit = sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])# 所有剩余數據集的誤差之和#print "sseSplit and notSplit: ",sseSplit,',',sseNotSplitif (sseSplit + sseNotSplit) < lowestSSE: # 劃分后的誤差和小于當前的誤差，本次劃分被保存#print 'here..........'bestCentToSplit = i # i代表簇數bestNewCents = centroidMat # 保存簇矩陣#print 'bestNewCents',bestNewCentsbestClustAss = splitClustAss.copy() # 拷貝所有數據點的簇索引和誤差lowestSSE = sseSplit + sseNotSplit # 保存當前誤差和# centList是原劃分的簇向量，bestCentToSplit是i值#print 'len(centList) and bestCentToSplit ',len(centList),',',bestCentToSplit# 數組過濾得到的是新劃分的簇類別是1的數據集的類別簇重新劃為新的類別值為最大的類別數bestClustAss[nonzero(bestClustAss[:,0].A == 1)[0],0] = len(centList) # 數組過濾得到的是新劃分的簇類別是0的數據集的類別簇重新劃為新的類別值為ibestClustAss[nonzero(bestClustAss[:,0].A == 0)[0],0] = bestCentToSplit#print 'the bestCentToSplit is: ',bestCentToSplit # 代表的是劃分的簇個數-1#print 'the len of bestClustAss is: ', len(bestClustAss) # 數據簇的數據點個數# 新劃分簇矩陣的第0簇向量新增到當前的簇列表中centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0] #print 'centList[bestCentToSplit]:',centList[bestCentToSplit]# 新劃分簇矩陣的第1簇向量添加到當前的簇列表中centList.append(bestNewCents[1,:].tolist()[0]) # centList是列表的格式#print 'centList',centList# 數組過濾得到所有數據集中簇類別是新簇的數據點clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]= bestClustAssreturn mat(centList), clusterAssment # 返回質心列表和簇分配結果# 球面距離計算，這里是利用球面余弦定理 def distSLC(vecA, vecB): # 經度和緯度用角度作為單位，這里用角度除以180然后乘以pi作為余弦函數的輸入a = sin(vecA[0,1]*pi/180) * sin(vecB[0,1]*pi/180) b = cos(vecA[0,1]*pi/180) * cos(vecB[0,1]*pi/180) * \cos(pi * (vecB[0,0]-vecA[0,0]) /180)return arccos(a + b)*6371.0 # 返回地球表面兩點之間的距離import matplotlib import matplotlib.pyplot as plt # 及簇繪圖函數 def clusterClubs(numClust=5): # 希望分得的簇數datList = [] # 創建一個空列表for line in open('places.txt').readlines():lineArr = line.split('\t')datList.append([float(lineArr[4]), float(lineArr[3])]) # 對應的是緯度和經度datMat = mat(datList) # 創建一個矩陣myCentroids, clustAssing = biKmeans(datMat, numClust, distMeas=distSLC)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') # imread()函數基于一幅圖像來創建矩陣ax0.imshow(imgP) # imshow()繪制該矩陣ax1=fig.add_axes(rect, label='ax1', frameon=False) # 在同一張圖上又創建一個字圖for i in range(numClust): # 遍歷每一個簇ptsInCurrCluster = datMat[nonzero(clustAssing[:,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(myCentroids[:,0].flatten().A[0], myCentroids[:,1].flatten().A[0], marker='+', s=300)plt.show() # 主函數 clusterClubs(5)

由于源代碼里給出了從Yahoo! placeFinder API得到的數據，所以這里直接拿來用了，得到的結果：

也可以通過改變聚類簇數來查看聚類效果，在此不再演示。

4. 注意的幾點：

4.1球面距離公式是計算球面上兩點間距離的公式

設所求
點A ：緯度β1 ，經度α1 ；
點B：緯度β2，經度α2；
則距離S=R?arccos[cosβ1cosβ2cos（α1?α2）+sinβ1sinβ2]
其中R<script type="math/tex" id="MathJax-Element-8">R</script>為球體半徑。

4.2 fig.add_axes()和Subplot（）用法

ax3 = fig.add_axes([0.1, 0.1, 0.8, 0.8]) ax4 = fig.add_axes([0.72, 0.72, 0.16, 0.16])

這里看幾個例子：

Axes - Subplot - Axis 之間的關系：

figure就是畫板，是畫紙的載體，但是具體畫畫等操作是在畫紙上完成的。在pyplot中，畫紙的概念對應的就是Axes/Subplot。

（1）subplot的操作

import matplotlib import matplotlib.pyplot as plt fig = plt.figure() ax = fig.add_subplot(111) ax.set(xlim=[0.5, 4.5], ylim=[-2, 8], title='An Example Axes', ylabel='Y-Axis', xlabel='X-Axis') plt.show()

運行結果：

（2）Axes 和 Subplot 的區別

Axes的使用：

fig = plt.figure() ax3 = fig.add_axes([0.1, 0.1, 0.8, 0.8]) ax4 = fig.add_axes([0.72, 0.72, 0.16, 0.16]) plt.show() print type(ax3)

運行結果：

這里軸域（Axes）可以理解成一些軸（Axis）的集合，當然這個集合還有很多軸（Axis）的屬性，標注等等。

我們用add_axes()方法生成一個軸域（Axes），括號里面的值前兩個是軸域原點坐標（這里指的是以整個figure為參照，前兩個值是相對于figure的左下角而言的，而不是當前子圖的起始坐標），后兩個是顯示坐標軸的長度（這里指的是子圖的軸長，也是相對于figure的長度而言的）。當我們生成了軸域的時候，從結果上看確實是生成了一個可以畫圖的子圖。

我們還可以分別在兩個軸域（Axes）中畫圖。 對比兩種方法，兩種對象，我們可以總結總結： add_subplot()方法在生成子圖過程，簡單明了，而用add_axes()方法,則生成子圖的靈活性更強，完全可以實現add_subplot()方法的功能，可以控制子圖顯示位置，甚至實現相互重疊的效果。下面是幾個例子：

subplot

fig = plt.figure() ax1 = fig.add_subplot(211) ax2 = fig.add_subplot(212) plt.show()

Axes

fig = plt.figure() ax3 = fig.add_axes([0, 0, 0.5, 0.5], title='An Example Axes', ylabel='Y-Axis', xlabel='X-Axis') ax4 = fig.add_axes([0.6, 0.6, 0.2, 0.16], title='An Example Axes', ylabel='Y-Axis', xlabel='X-Axis') plt.show()

可以看出ax3是從畫板的原點（0，0）開始畫，橫縱坐標軸都是0.5
ax4是從畫板的（0.6，0.6）開始畫，橫縱坐標軸分別是0.2和0.6

fig = plt.figure() ax3 = fig.add_axes([0, 0, 0.3, 1], title='An Example Axes', ylabel='Y-Axis', xlabel='X-Axis') ax4 = fig.add_axes([0.6, 0.3, 0.2, 0.3], title='An Example Axes', ylabel='Y-Axis', xlabel='X-Axis') plt.show()

參考：https://www.zhihu.com/question/51745620

4.3 scatterMarkers[i % len(scatterMarkers)]

markerStyle = scatterMarkers[i % len(scatterMarkers)]

使用：

>>> scatterMarkers=['s', 'o', '^', '8', 'p', 'd', 'v', 'h', '>', '<'] >>> for i in range(5): ... markerStyle = scatterMarkers[i % len(scatterMarkers)] ... print markerStyle ... s o ^ 8 p >>>

類似于列表推到式，但此還有一個特殊的地方是可以循環使用，也就是更多簇時，可以循環標記。

總結

以上是生活随笔為你收集整理的K-均值对地图上的点进行聚类（2）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。