Python地信专题 |基于geopandas的空间数据分析-深入浅出分层设色
點擊藍字關注我，有干貨領取！
本文對應代碼和數據已上傳至我的Github倉庫：
https://github.com/CNFeffery/DataScienceStudyNotes[1]
已發布：
Python地信專題 | 基于geopandas的空間數據分析—數據結構篇
Python地信專題 | 基于geopandas的空間數據分析-坐標參考系篇
Python地信專題 | 基于geopandas的空間數據分析-文件IO篇
Python地信專題 | 基于geopandas玩轉地圖可視化
1 簡介
通過前面的文章，我們已經對geopandas中的數據結構、坐標參考系、文件IO以及基礎可視化有了較為深入的學習。
其中在基礎可視化那篇文章中我們提到了分層設色地圖，可以對與多邊形關聯的數值屬性進行分層，并分別映射不同的填充顏色。
但只是開了個頭舉了個簡單的例子，實際數據可視化過程中的分層設色有一套策略方法。
作為基于geopandas的空間數據分析系列文章的第五篇，通過本文你將會學習到基于geopandas和機器學習的分層設色。
2 基于geopandas的分層設色
地區分布圖（Choropleth maps，又叫面量圖）作為可能是最常見的一種地理可視化方法。
其核心是對某個與矢量面關聯的數值序列進行有意義的分層，并為這些分層選擇合適美觀的色彩，最后完成對地圖的著色。
優點是美觀且直觀，即使對地理信息一竅不通的人，也能通過顏色區分出不同面之間的同質性與異質性：
圖1但同樣地，如果對數據分層采取的方法有失嚴謹沒有很好的遵循數據特點，會很容易讓看到圖的人產生出不正確的判斷。
下面我們按照先分層，后設色的順序進行介紹。
2.1 基于mapclassify的數據分層
上一篇文章中我們提到過,，在geopandas.GeoDataFrame.plot()中，參數scheme對應的數據分層是基于第三方庫mapclassify實現的。
因此要想對geopandas中的數據分層有深入的了解，我們就得先來了解一下mapclassify中的各種數據分層算法。
用到的數據是系列文章前幾期使用地滾瓜爛熟的新冠肺炎疫情數據，數據處理過程同上一篇文章，這里不再解釋：
圖22.1.1 BoxPlot
BoxPlot即箱線圖，是統計學中使用到的一種方法：
對個數為? 觀測數據從小到大進行排序，分別得到位置處于?、 以及? 的觀測值，稱為?、 以及?（即第一四位數、中位數和第三四分位數）。
并定義? 為?，以? 為下限，以? 為上限，將小于下限或大于上限的觀測值作為離群異常值。
最后用圖像的形式表達上述計算結果，如圖2的上圖，而圖2的下圖對應著概率估計。
可以看出，箱線圖法實際上是基于概率估計的一種異常值剔除方法，因為離群值只有? 的概率會出現，即如果你想要找出數據中的異常高低值，BoxPlot是不錯的選擇：
圖3在mapclassify中我們使用BoxPlot()來為數據實現箱線圖分層：
import mapclassify as mc# 對各省2020-03-08對應的累計確診數量進行分層 bp = mc.BoxPlot(temp['province_confirmedCount']) # 查看數據分層結果 bp 圖4可以看出通過箱線圖法將數據分成了五類，其中異常值只有1個即為湖北省。
下面我們配合geopandas來對上述結果進行可視化，和上一篇文章一樣，按照省級單位名稱連接我們的疫情數據與矢量數據：
圖5接著對其進行可視化，在上一篇文章圖28的基礎上，將scheme參數改為BoxPlot，又因為箱線圖可以看作無監督問題，故分層數量k在這里無效，刪去：
fig, ax = plt.subplots(figsize=(10, 10))ax = data_with_geometry.to_crs(albers_proj).plot(ax=ax,column='province_confirmedCount',cmap='Reds',missing_kwds={"color": "lightgrey","edgecolor": "black","hatch": "","label": "缺失值"},legend=True,scheme='BoxPlot',legend_kwds={'loc': 'lower left','title': '確診數量分級','shadow': True})ax = nine_lines.geometry.to_crs(albers_proj).plot(ax=ax,edgecolor='grey',linewidth=3,alpha=0.4)ax.axis('off') plt.suptitle('新型冠狀肺炎累計確診數量地區分布', fontsize=24) # 添加最高級別標題 plt.tight_layout(pad=4.5) # 調整不同標題之間間距 ax.text(-2800000, 1300000, '* 原始數據來源：丁香園，\n其中臺灣及香港數據缺失') # 添加數據說明fig.savefig('圖6.png', dpi=300) 圖6咋看起來沒問題，但是如果你仔細觀察左下角的圖例會發現前兩行范圍顏色是重復的，且數值范圍是錯亂的。
這是geopandas.GeoDataFrame.plot()中涉及箱線圖法的一個小bug，遇到這種問題不用慌。
如果你在上一篇文章中去我的Github倉庫查看過創作圖29對應的代碼，一定會想到既然geopandas自身有bug，那我們用matplotlib中的mpatches和legend自定義圖例就可以啦。
而為了自定義的圖例色彩與geopandas映射出的保持一致，我們需要額外使用到matplotlib中的get_cmap(cmap)來制作可獨立導出顏色的cmap方案實例。
譬如我們這里是Reds，就需要按照前面bp的有記錄數量的分層結果，從Reds中產生同樣5個檔次的顏色，具體操作過程如下：
import matplotlib.patches as mpatchesfig, ax = plt.subplots(figsize=(10, 10))ax = data_with_geometry.to_crs(albers_proj).plot(ax=ax,column='province_confirmedCount',cmap='Reds',missing_kwds={"color": "lightgrey","edgecolor": "black","hatch": "","label": "缺失值"},scheme='BoxPlot')handles, labels = ax.get_legend_handles_labels() #get existing legend item handles and labelsax = nine_lines.geometry.to_crs(albers_proj).plot(ax=ax,edgecolor='grey',linewidth=3,alpha=0.4)# 實例化cmap方案 cmap = plt.get_cmap('Reds')# 得到mapclassify中BoxPlot的數據分層點 bp = mc.BoxPlot(temp['province_confirmedCount']) bins = bp.bins# 制作圖例映射對象列表 LegendElement = [mpatches.Patch(facecolor=cmap(_*0.25), label=f'{int(max(bins[_], 0))} - {int(bins[_+1])}')for _ in range(5)] + \[mpatches.Patch(facecolor='lightgrey', edgecolor='black', hatch='', label='缺失值')]# 將制作好的圖例映射對象列表導入legend()中，并配置相關參數 ax.legend(handles = LegendElement, loc='lower left', fontsize=10, title='確診數量分級', shadow=True, borderpad=0.6)ax.axis('off') plt.suptitle('新型冠狀肺炎累計確診數量地區分布(截至2020年03月04日)', fontsize=24) # 添加最高級別標題 plt.title('數據分層方法：BoxPlot', fontsize=18) plt.tight_layout(pad=4.5) # 調整不同標題之間間距 ax.text(-2900000, 1250000, '* 原始數據來源：丁香園，\n其中臺灣及香港數據缺失') # 添加數據說明fig.savefig('圖7.png', dpi=300) 圖7可以看到，通過自定義圖例的方式，雖然麻煩了一點，但是我們不僅修復了圖例的bug，還為其添加了更加完善的細節。
如圖形修改為矩形，范圍修改為整數。
2.1.2 EqualInterval
EqualInterval即等間距，是最簡單的一種分層方法。
它在原數據最小值與最大值間以等間距的方式劃分出k個層次，mapclassify中對應等間距法的類為EqualInterval()：
bp = mc.EqualInterval(temp['province_confirmedCount']) # 查看數據分層結果 bp 圖8可以看到對于分布非常不均勻的新冠肺炎確診數量數據來說，這種方法表現得十分糟糕，中間三個類都沒有記錄落入。
如果使用這種方法強行繪圖，效果就會類似上一篇文章中地區分布圖部分。
最開始那個糟糕的效果那樣只有湖北一個地方是最深的暗紅色，而其他地方皆為最淡的色階，這里就不重復演示。
2.1.3 FisherJenks
在了解mapclassify中的FisherJenks之前，我們先來了解一下什么是Jenks Natural Breaks：
Jenks Natural Breaks
Jenks Natural Breaks旨在為1維數據計算合適的劃分點，使得不同組之間的差距盡可能大的同時組內差距盡可能小。
其思路非常簡單，舉一個簡單的例子進行說明：
對于一組待分割的序列?，現在需要為其找到將原始數據分為? 部分的方法。
那么實際上就有?、 以及? 這三種切分方法，現定義sum of squared deviations for array mean（簡稱SDAM）：
以及針對每一種數據分層方法，在其分出的每一組? 上計算組內離差平方和并累加所有組的結果，定義為sum of squared deviations for class means（簡稱SDCM_ALL）：
有了? 和?，現在對分組優劣定義一個評判指標goodness of variance fit（簡稱GVF），取值范圍為?，越高越好：
這樣我們就可以對每一種分組方案進行評價，譬如對我們上面簡單的例子：
則對應各種方案的GVF計算如下：
可以看出第二種方案? 的分層方法效果最好，也與我們對數據的直觀感覺相貼合。
這就是Jenks Natural Breaks的基本思路，但這種暴力遍歷所有分組方案的做法對數據數量及選擇分組的個數很敏感。
尤其是對分組數量，一旦分組數量過于多，待篩選計算的方案數量就變成了天文數字。
下面我來告訴大家為什么：
定義長度為? 的序列?。
且滿足? 時?，即整個序列從小到大單調遞增，那么將其分成? 組的過程，可以分解為先選擇第一組。
且為了保證右邊剩余? 個組每組至少有1個數據分配，則第一組有? 種分配方式，而第一組包含的數字數量? 確定之后，剩余? 個數據的繼續分組又可以視為獨立的遞歸分組過程。
因此最終需要考慮的方案個數用公式表達起來有些復雜，但是換成計算機中的遞歸過程就變得一目了然。
我經過思考和紙上的推演，寫出了下面所示的遞歸函數f(n, k)來實現方案總數的計算：
def f(n, k):# 若k退化為2，則顯然需要n - 1種方案，譬如4個數字分2組有3種方案if k == 2:return n - 1else:# 若k未退化為2，則繼續遞歸過程return sum([f(n-_, k-1) for _ in range(1, n - k + 2)])有了這個遞歸函數，我們就可以來直觀的看一看為什么不能選擇太多分組。
首先我們對長度為100的序列分為3組試試：
f(100, 3) Out[11]: 4851可以看到待選擇的方案才4851個，還是很少的，那么我們接下來將組數提高到5：
f(100, 5) Out[12]: 3764376發生了什么？隨著遞歸深度的增大，待選擇方案數量一下子就提高到三百多萬個！再切換成7試一下：
f(100, 7) Out[13]: 1120529256在跑上述代碼時，明顯能感受到計算花費時間的激增，最終結果也達到驚人的11億多！
看到這，我們就明白了，原始的Jenks Natural Breaks算法雖然很有效，但如果以暴力遍歷的方式計算，其復雜度是難以應付日常需求的。
為了對其進行優化，以在少量的計算時間內計算出盡可能靠譜的分組結果，一系列改良加速方法被提出，而mapclassify中的FisherJenks，即為jenks教授在論文Fisher, W. D., 1958, On grouping for maximum homogeneity.的基礎上提出的改良算法。
但這是一個很神秘的算法，根據Tom MacWright文章[2] 介紹，jenks教授的原始論文沒有留下數字化資料，一直為堪薩斯大學地理學系所私有。
而隨著1996年jenks教授的離世，原論文需要到2072年版權才能到期公開，所以我們現在在各種GIS類軟件以及各種開源軟件包中使用到的fisher jenks算法，均是對最初的一段Fortran代碼的移植和改造，這也成了一段未解之謎。
感興趣的讀者可以去文章Natural Breaks Classification[3] 了解更多。
回到我們的主題，搞清楚了FisherJenks的計算目標之后，我們同樣利用mapclassify計算分層結果，其默認分層為5：
圖9可以看到，在這種方式下，數據的分組較為合理，同樣將geopandas.GeoDataFrame.plot()中的參數設置為FisherJenks繪制出圖10：
圖10與BoxPlot相比差距還是比較明顯，處于第二級嚴重程度的省份只有河南、廣東及浙江，更貼近數據的自然層次結構。
2.1.4 NaturalBreaks
等下！上一小結中的FisherJenks不就是我們俗稱的自然斷點法嗎，怎么又來了個NaturalBreaks？
其實我在翻看mapclassify的官方文檔看到這里時，也很疑惑。
于是我仔細研究了NaturalBreaks對應的源代碼，追根溯源，WHAT?，竟然是k-mean算法，而且直接調用的scikit-learn的KMEANS。。。
圖11不過也可以理解，畢竟k-means就是在找數據中組內相似度盡可能高且組間差異盡量大的簇，關于k-means我想我就不需要贅述了，畢竟是最基礎的數據挖掘算法之一。
而scikit-learn里默認的KMEANS使用的k-means++初始方式，只是在原始k-means基礎上，修改了后續初始點的概率密度，使得k-means算法更加魯棒穩定。
下面直接來看NaturalBreaks的數據分層結果：
圖12和FisherJenks的結果竟然一樣，但如果你多運行幾次會發現這個結果不是完全固定的，由于k-means隨機初始迭代起點。
因此不同次運行的結果可能會有輕微差別（圖13），在數據量很大時，基于快速聚類法的NaturalBreaks是較為理想的數據分層選擇：
圖13配合geopandas繪圖只需要把scheme參數修改為NaturalBreaks即可，因為跟FisherJenks類似，這里就不再贅述。
2.1.5 JenksCaspall
mapclassify中的JenksCaspall本質上為k-medians聚類，其首先根據分層層數? 在數據中找到? 個分位數點，將原始數據等分為數量盡可能相同的? 份并以這? 份數據的中位數作為各自的初始點。
接著基于k-medians的思想，迭代計算為每個樣本點找到與其距離更近的中位數點，并以此重新劃分分層以及重新計算各分層中位數點，直至每個數據對應的分層標簽不再變化。
再將每個分層中數據的最大值作為間斷點，下面我們從mapclassify源代碼中抽出該部分代碼，對其迭代過程可視化。
具體的代碼較多，請在文章開頭的Github倉庫中對應本文路徑下查看：
圖14其中顏色區分對應迭代輪次的數據分層歸屬，虛線代表對應迭代輪次的間斷點，仔細可以看出在迭代過程中數據分層的變化情況。
用JenksCaspall數據分層出來的結果，無論數據分布如何，每個分層內部的數據個數都較為均勻，下面我們用JenksCaspall來劃分省份疫情嚴重情況：
圖15可以看到被分到最嚴重級別的不再只有湖北省，當你希望數據分層個數較為均勻時，JenksCaspall是個不錯的選擇。
2.1.6 HeadTailBreaks
HeadTailBreaks是一種較為嶄新的數據分層方法，出自*Head/Tail Breaks: A New Classification Scheme for Data with a Heavy-Tailed Distribution*[4]，專門用于對具有重尾特點的數據進行分層。
所謂重尾即在整個數據中，較小的值數量往往較多，而最大的位于頭部的值數量很少，其數據分布呈現出“尾重頭輕”的特點：
圖16這種典型如人口密度分布數據，數值較低的點往往數量眾多，聚集在尾部，形成重尾，HeadTailBreaks的優點是可以盡量在地區分布圖中真實反映原始數據的分布特點。
如圖17（https://sites.google.com/site/thepowerofcartography/head-tail-breaks），左邊是FisherJenks，右邊是HeadTailBreaks，可以看出，右圖相對于左圖更好地體現了原始數據的重尾特點，最淺色的圖斑數量明顯多于次淺色的圖斑：
圖17在geopandas中使用時傳入scheme='HeadTailBreaks'即可（由于新冠肺炎各省份確診數量數據尾部和頭部最大值之間沒有較為連續的中間值過渡，不太適合用此方法故不作演示）。
2.1.7 Quantiles
Quantiles即分位數，原理很簡單，根據分位數點對原數據進行等分：
圖18利用Quantiles對確診數量分組可視化：
圖192.1.8 Percentiles
同樣是使用分位數對數據進行分層，Percentiles提供了參數pct以允許用戶以百分位數的形式傳入自定義分隔點。
譬如我們將[1, 50, 99, 100]作為pct的傳入值，則分組結果如下：
圖20每個傳入的百分位點其左邊到上一個分隔點為止，包括其本身，將被分到同一組。
對應的圖像如圖21，在geopandas中使用時除了設置scheme='Percentiles'之外，還要在另一個字典型參數classification_kwds中傳入{'pct': 百分位數列表}：
圖212.1.9 StdMean
StdMean的思想類似前面的箱線圖，不同的是箱線圖屬于非參數方法，而StdMean建立在正態分布為基礎的經驗法則之上。
即對于正態分布而言，68%的數據將分布在距離均值1個標準差之內，95%的數據在2個標準差之內，99.7%的數據在3個標準差之內。
即對原始數據標準化之后，根據距離樣本均值的不同標準差范圍來劃分數據，mapclassify中的StdMean默認按照[-2, -1, 1, 2]來劃分：
圖222.1.10 UserDefined
關于數據分層最后要介紹的是自定義分層，即按照用戶輸入的分隔點來自由劃分數據集，譬如我們按照新浪新聞疫情地圖的劃分方式：
圖23 結合`geopandas`使用時除了設置`scheme='UserDefined'`以外，還要設置`classification_kwds`中的`bins=分隔點列表`：fig, ax = plt.subplots(figsize=(10, 10))ax = data_with_geometry.to_crs(albers_proj).plot(ax=ax,column='province_confirmedCount',cmap='Reds',missing_kwds={"color": "lightgrey","edgecolor": "black","hatch": "","label": "缺失值"},legend=True,scheme='UserDefined',classification_kwds={'bins': [9, 99, 499, 999, 9999]},legend_kwds={'loc': 'lower left','title': '確診數量分級','shadow': True})ax = nine_lines.geometry.to_crs(albers_proj).plot(ax=ax,edgecolor='grey',linewidth=3,alpha=0.4)ax.axis('off') plt.suptitle('新型冠狀肺炎累計確診數量地區分布', fontsize=24) # 添加最高級別標題 plt.tight_layout(pad=4.5) # 調整不同標題之間間距 ax.text(-2800000, 1300000, '* 原始數據來源：丁香園，\n其中臺灣及香港數據缺失') # 添加數據說明fig.savefig('圖24.png', dpi=300) 圖242.2 色彩方案的選擇
前面已經詳細介紹了數據分層常用的各種方法及使用場景，“分層”的部分做完之后，就到了設色的部分。
其實色彩搭配是比較主觀的事情，但想要自己創造出美觀合理的配色方案并不是容易的事情。
下面我們來介紹兩種選擇配色方案的方法。
2.2.1 基于palettable的配色
下面我要給大家介紹的Python第三方庫palettable在我之前關于詞云圖的一篇文章中介紹stylecloud時介紹過，是專門幫助我們為可視化作品配色的。
palettable不依賴其他三方庫，純Python實現，其強大之處在于內置了數量驚人的經典配色方案[5]。
囊括了CartoColors、cmocean、Colorbrewer2、Cubehelix、Light & Bartlein、matplotlib、MyCarta、Scientific、Tableau以及The Wes Anderson Palettes blog中的大量經典配色方案：
-palettable.cartocolors.diverging
-palettable.cartocolors.qualitative
-palettable.cartocolors.sequential
-palettable.cmocean.diverging
-palettable.cmocean.sequential
-palettable.colorbrewer.diverging
-palettable.colorbrewer.qualitative
-palettable.colorbrewer.sequential
-palettable.lightbartlein.diverging
-palettable.lightbartlein.sequential
-palettable.matplotlib
-palettable.mycarta
-palettable.scientific.diverging
-palettable.scientific.sequential
-palettable.tableau
-palettable.wesanderson
使用起來非常簡單。
譬如如果我們想要使用palettable.cmocean.sequential中的色彩，其中cmocean表示色彩來源，sequential表示連續型色彩，就可以先在對應的示例網頁[6]下查看所有方案：
圖25比如我對其中的Dense方案很中意：
圖26就可以按照如下方式，先從palettable中導入對應顏色。
譬如我們導入Dense_20，20表示其自帶的離散色彩數量，并查看其自帶的離散色彩RGB值、離散色盤以及連續色盤示例：
from palettable.cmocean.sequential import Dense_20 from pprint import pprintprint('對應離散顏色：') pprint(Dense_20.colors) print('離散：') Dense_20.show_discrete_image() print('連續:') Dense_20.show_continuous_image() 圖27使用.mpl_colormap將其轉換為matplotlib可接受的cmap數據結構，作為cmap參數值傳入繪圖部分即可：
圖28如果想要翻轉映射方向，換成Dense_20_r再重復上述操作即可：
圖29更多palettable自帶色彩方案，可以在https://jiffyclub.github.io/palettable/ 下查看探索。
2.2.2 基于圖片主色的配色
我們在生活中偶然會看到配色方案讓人眼前一亮的海報或畫作，這時如果你想將這些作品中的主要顏色也應用到自己的可視化作品上，可以參考我下面的做法。
這里以我很喜歡的賈樟柯導演的《一直游到海水變藍》中文版海報為例：
圖30思路是抽取所有像素點的RGB三通道值，分別作為三個特征，輸入k-means中進行聚類，將聚類數量設置為你想要提取出的主色數量：
from sklearn.cluster import KMeans# 構建特征 rgb = pd.DataFrame([sea[x][y] for x in range(sea.shape[0]) for y in range(sea.shape[1])],columns=['r', 'g', 'b'])# k-means聚類，其中n_clusters表示聚類數量，n_jobs=-1表示開啟所有核心并行運算 model = KMeans(n_clusters=5, n_jobs=-1) model.fit(rgb) # 訓練模型# 提取聚類簇重心，即我們需要的主色，繪制調色板 plt.bar([i for i in range(model.cluster_centers_.__len__())],height=[1 for i in range(model.cluster_centers_.__len__())],color=[tuple(c) for c in (model.cluster_centers_ / 255.)],width=1) plt.axis('off') 圖31再來個例子，提取《一直游到海水變藍》海外版海報主色：
圖32對應提取到的5種主色如圖33：
圖33類似的，你可以試著提取你喜愛的平面作品的主色。
以上就是本文的全部內容，如有筆誤望指出。
參考資料
[1]
https://github.com/CNFeffery/DataScienceStudyNotes: https://github.com/CNFeffery/DataScienceStudyNotes
[2]Tom MacWright文章: https://macwright.org/2013/02/18/literate-jenks.html
[3]Natural Breaks Classification: https://stat.ethz.ch/pipermail/r-sig-geo/2006-March/000811.html
[4]Head/Tail Breaks: A New Classification Scheme for Data with a Heavy-Tailed Distribution: https://www.tandfonline.com/doi/abs/10.1080/00330124.2012.700499
[5]經典配色方案: https://jiffyclub.github.io/palettable/cartocolors/diverging/
[6]示例網頁: https://jiffyclub.github.io/palettable/cmocean/sequential/
-END-
往期精彩回顧適合初學者入門人工智能的路線及資料下載機器學習在線手冊深度學習在線手冊AI基礎下載（pdf更新到25集）本站qq群1003271085，加入微信群請回復“加群”獲取一折本站知識星球優惠券，請回復“知識星球”喜歡文章，點個在看
總結
以上是生活随笔為你收集整理的Python地信专题 |基于geopandas的空间数据分析-深入浅出分层设色的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 肖战被抵制？Python爬虫揭秘关于肖战
- 下一篇： 图解Reformer：一种高效的Tran