（資料科學學習手札83）基於geopandas的空間資料分析——geoplot篇(下)
阿新 • • 發佈：2020-05-16
> 本文示例程式碼、資料及檔案已上傳至我的`Github`倉庫[https://github.com/CNFeffery/DataScienceStudyNotes](https://github.com/CNFeffery/DataScienceStudyNotes)
# 1 簡介
在上一篇文章中我們詳細學習了`geoplot`中較為基礎的三種繪圖API：`pointplot()`、`polyplot()`以及`webmap()`，而本文將會承接上文的內容，對`geoplot`中較為實用的幾種高階繪圖API進行介紹。
圖1
本文是***基於geopandas的空間資料分析***系列文章的第7篇，通過本文你將學習`geoplot`中的高階繪圖API。
# 2 geoplot進階
上一篇文章中的`pointplot()`、`polyplot`以及`webmap()`幫助我們解決了在繪製散點、基礎面以及新增線上地圖底圖的問題，為了製作出資訊量更豐富的視覺化作品，我們需要更強的操縱向量資料與對映值的能力，`geoplot`為我們封裝好了幾種常見的高階視覺化API。
## 2.1 Choropleth
*Choropleth*圖又稱作**地區分佈圖**或**面量圖**，我們在系列之前的***深入淺出分層設色篇***中介紹過其原理及`geopandas`實現，可以通過將指標值對映到面數據上，以實現對指標值地區分佈的視覺化。
在`geoplot`中我們可以通過`choropleth()`來快速繪製地區分佈圖，其主要引數如下：
> 
圖2
圖3
首先我們將兩張表中各自對應的州名資料作為鍵進行連線（注意`pd.merge`返回的結果型別為`DataFrame`，需要轉換回`GeoDataFrame`）：
```Python
# 按照州名列進行連線
usa_plot_base = pd.merge(left=contiguous_usa,
right=usa_covid19_20200513,
left_on='state',
right_on='Province_State')
# 轉換DataFrame到GeoDataFrame，注意要加上crs資訊
usa_plot_base = gpd.GeoDataFrame(usa_plot_base, crs='EPSG:4326')
```
接下來我們將確診數作為對映值，因為美國各州中紐約州和新澤西州確診數量分別達到了34萬和14萬，遠遠超過其他州，所以這裡作為單獨的圖層進行陰影填充以突出其嚴重程度：
```Python
# 圖層1：除最嚴重兩州之外的其他州
ax = gplt.choropleth(df=usa_plot_base.query("state not in ['New York', 'New Jersey']"),
projection=gcrs.AlbersEqualArea(),
hue='Confirmed',
scheme=mc.FisherJenks(usa_plot_base.query("state not in ['New York', 'New Jersey']")['Confirmed'], k=3),
cmap='Reds',
alpha=0.8,
edgecolor='lightgrey',
linewidth=0.2,
figsize=(8, 8)
)
# 圖層2：紐約州
ax = gplt.polyplot(df=usa_plot_base.query("state == 'New York'"),
hatch='/////',
edgecolor='black',
ax=ax)
# 圖層3：新澤西州
ax = gplt.polyplot(df=usa_plot_base.query("state == 'New Jersey'"),
hatch='/////',
edgecolor='red',
extent=usa_plot_base.total_bounds,
ax=ax)
# 例項化cmap方案
cmap = plt.get_cmap('Reds')
# 得到mapclassify中BoxPlot的資料分層點
bp = mc.FisherJenks(usa_plot_base.query("state not in ['New York', 'New Jersey']")['Confirmed'], k=3)
bins = [0] + bp.bins.tolist()
# 製作圖例對映物件列表，這裡分配Greys方案到三種色彩時對應的是[0, 0.5, 1]這三個取樣點
LegendElement = [mpatches.Patch(facecolor=cmap(_ / 2), label=f'{int(bins[_])}-{int(bins[_+1])}')
for _ in range(3)] + \
[mpatches.Patch(facecolor='none',
edgecolor='black',
linewidth=0.2,
hatch='/////',
label='New York: {}'.format(usa_plot_base.query("state == \"New York\"").Confirmed.to_list()[0])),
mpatches.Patch(facecolor='none',
edgecolor='red',
linewidth=0.2,
hatch='/////',
label='New Jersey: {}'.format(usa_plot_base.query("state == \"New Jersey\"").Confirmed.to_list()[0]))]
# 將製作好的圖例對映物件列表匯入legend()中，並配置相關引數
ax.legend(handles = LegendElement,
loc='lower left',
fontsize=8,
title='確診數量',
title_fontsize=10,
borderpad=0.6)
# 新增標題
plt.title('美國新冠肺炎各州病例數(截至2020.05.14)', fontsize=18)
# 儲存影象
plt.savefig('圖4.png', dpi=300, pad_inches=0, bbox_inches='tight')
```
圖4
這樣我們就得到了圖4，需要注意的是，`geoplot.choropleth()`只能繪製地區分佈圖，傳入面數據後`hue`引數必須指定對應對映列，否則會報錯，因此這裡我們疊加紐約州和新澤西州單獨面圖層時使用的是`polyplot()`。
## 2.2 Kdeplot
`geoplot`中的`kdeplot()`對應核密度圖，其基於`seaborn`中的`kdeplot()`，通過對向量點資料分佈計算核密度估計，從而對點資料進行視覺化，可用來展示點資料的空間分佈情況，其主要引數如下：
> **df**：傳入對應的存放點資料的`GeoDataFrame`物件
>
> **projection**：用於指定投影座標系，傳入`geoplot.crs`中的物件
>
> **cmap**：和`matplotlib`中的**cmap**使用方式一致，用於控制色彩對映方案
>
> **clip**：`GeoSeries`型，用於為初始生成的核密度影象進行蒙版裁切，下文會舉例說明
>
> **extent**：元組型，用於傳入左下角和右上角經緯度資訊來設定地圖空間範圍，格式為`(min_longitude, min_latitude, max_longitude, max_latitude)`
>
> **figsize**：元組型，用於控制畫幅大小，格式為`(x, y)`
>
> **ax**：`matplotlib`座標軸物件，如果需要在同一個座標軸內疊加多個圖層就需要用這個引數傳入先前待疊加的`ax`
>
> **shade**：bool型，當設定為False時只有等值線被繪製出，當設定為True時會繪製核密度填充
>
> **shade_lowest**：bool型，控制是否對概率密度最低的層次進行填充，下文會舉例說明
>
> **n_levels**：int型，控制等值線數量，即按照概率密度對空間進行均勻劃分的數量
下面我們回到上一篇文章開頭的例子——紐約車禍記錄資料，在其他引數均為預設的情況下，呼叫`kdeplot`對車禍記錄點資料的空間分佈進行視覺化：
```Python
# 圖層1：行政邊界
ax = gplt.polyplot(df=nyc_boroughs,
projection=gcrs.AlbersEqualArea())
# 圖層2：預設引數下的kdeplot
ax = gplt.kdeplot(df=nyc_collision_factors,
cmap='Reds',
ax=ax)
# 儲存影象
plt.savefig('圖5.png', dpi=300, pad_inches=0, bbox_inches='tight')
```
圖5
可以看到，在`kdeplot()`主要引數均為預設值的情況下，我們得到了點資料空間分佈的概率估計結果及其等高線，譬如圖中比較明顯能看到的兩個點分佈較為密集的中心，下面我們調整`n_levles`引數到比較大的數字：
```Python
# 圖層1：行政邊界
ax = gplt.polyplot(df=nyc_boroughs,
projection=gcrs.AlbersEqualArea())
# 圖層2：kdeplot
ax = gplt.kdeplot(df=nyc_collision_factors,
cmap='Reds',
n_levels=30,
ax=ax,
figsize=(8, 8))
# 儲存影象
plt.savefig('圖6.png', dpi=300, pad_inches=0, bbox_inches='tight')
```
圖6
可以看到在增大`n_levels`引數後，圖中等值線的數量隨之增加，下面我們設定`shade=True`：
```Python
# 圖層1：行政邊界
ax = gplt.polyplot(df=nyc_boroughs,
projection=gcrs.AlbersEqualArea())
# 圖層2：kdeplot
ax = gplt.kdeplot(df=nyc_collision_factors,
cmap='Reds',
shade=True,
ax=ax,
figsize=(8, 8))
# 儲存影象
plt.savefig('圖7.png', dpi=300, pad_inches=0, bbox_inches='tight')
```
圖7
這時影象等值線間得到相應顏色的填充，使得點分佈中心看起來更加明顯，再新增引數`shade_lowest=True`，即可對空白區域進行填充：
圖8
隨之而來的問題是整幅影象都被填充，為了裁切出核密度影象的地區輪廓，將底層行政區面數據作為`clip`的引數傳入，便得到理想的效果：
圖9
## 2.3 Sankey
桑基圖專門用於表現不同物件之間某個指標量的流動情況，譬如最常見的航線流向情況，其本質是對線資料進行視覺化，並將指標值對映到線的色彩或粗細水平上，而`geoplot`中的`sankey()`可以用來繪製這種圖，尷尬的是`sankey()`繪製出的OD流向圖實在太醜，但`sankey()`中將數值對映到線資料色彩和粗細的特性可以用來進行與流量相關的視覺化，其主要引數如下：
> **df**：傳入對應的`GeoDataFrame`物件
>
> **projection**：用於指定投影座標系，傳入`geoplot.crs`中的物件
>
> **hue**：傳入對應**df**中指定列名或外部序列資料，用於對映線的顏色，預設為None即不進行設色
>
> **cmap**：和`matplotlib`中的**cmap**使用方式一致，用於控制色彩對映方案
>
> **alpha**：控制全域性色彩透明度
>
> **scheme**：作用類似`geopandas`中的**scheme**引數，用於控制分層設色，詳見本系列文章的***分層設色篇***，但不同的是在`geoplot`0.4.0版本之後此引數不再搭配分層數量*k*共同使用，而是更新為傳入`mapclassify`分段結果物件，下文中會做具體演示
>
> **scale**：用於設定對映線要素粗細程度的序列資料，格式同**hue**，預設為None即每條線等粗
>
> **linewidth**：當不對線寬進行對映時，該引數用於控制線寬
>
> **legend**：bool型，用於控制是否顯示圖例
>
> **legend_values**：list型，用於自定義圖例顯示的各個具體數值
>
> **legend_labels**：list型，用於自定義圖例顯示的各個具體數值對應的文字標籤，與**legend_values**搭配使用
>
> **legend_kwargs**：字典，在**legend**引數設定為True時來傳入更多微調圖例屬性的引數
>
> **extent**：元組型，用於傳入左下角和右上角經緯度資訊來設定地圖空間範圍，格式為`(min_longitude, min_latitude, max_longitude, max_latitude)`
>
> **figsize**：元組型，用於控制畫幅大小，格式為`(x, y)`
>
> **ax**：`matplotlib`座標軸物件，如果需要在同一個座標軸內疊加多個圖層就需要用這個引數傳入先前待疊加的`ax`
下面我們以2015年華盛頓街道路網日平均交通流量資料為例，其中每個要素均為線要素，`aadt`代表日均流量：
圖10
我們將其流量列對映到線的粗細程度和顏色上來，為了美觀起見我們選擇系列文章**分層設色篇**中`palettable`的`SunsetDark`作為配色方案：
```Python
# 選擇配色方案為SunsetDark_5
from palettable.cartocolors.sequential import SunsetDark_5
gplt.sankey(
dc_roads,
projection=gcrs.AlbersEqualArea(),
scale='aadt',
hue='aadt',
limits=(0.1, 2), # 控制線寬範圍
scheme=mc.NaturalBreaks(dc_roads['aadt']),
cmap=SunsetDark_5.mpl_colormap,
figsize=(8, 8),
extent=dc_roads.total_bounds
)
plt.savefig("圖11.png", dpi=500, pad_inches=0, bbox_inches='tight')
```
圖11
## 2.4 geoplot中的座標參考系
`geoplot`中的座標參考系與`geopandas`中管理起來的方式截然不同，因為`geopandas`基於`pyproj`管理座標參考系，而`geoplot`中的`crs`子模組來源於`cartopy`，這一點我跟`geoplot`的主要開發者聊過，他表示`geoplot`暫時不支援`geopandas`中那樣自定義任意投影或使用`EPSG`投影，而是內建了一系列常用的投影，譬如我們上文中繪製美國區域時頻繁使用到的`AlbersEqualArea()`即之前我們在`geopandas`中通過proj4自定義的阿爾伯斯等面積投影，其他常見投影譬如*Web Mercator*、*Robinson*，或者直接繪製球體地圖，如本文開頭的圖1就來自官方示例（https://residentmario.github.io/geoplot/gallery/plot_los_angeles_flights.html#sphx-glr-gallery-plot-los-angeles-flights-py），關於`geoplot`座標參考系的細節比較簡單本文不多贅述，感興趣的讀者可以前往官網（https://residentmario.github.io/geoplot/api_reference.html#projections）檢視。
## 2.5 在模仿中學習
又到了最喜歡的“復刻”環節啦，本文要模仿的地圖視覺化作品來自https://github.com/Z3tt/30DayMapChallenge/tree/master/contributions/Day26_Hydrology，同樣是用`R`語言實現，對全球主要河流的形態進行優雅地視覺化：
圖12
針對其河流寬度方面的視覺化，我們基於上文中的`sankey()`來實現，由於原圖中南極洲區域實際上是誇大了的，其`R`原始碼中設定的緯度範圍達到了-110度，這是原作者為了放得下標題內容，所以在影象下部區域虛構了一篇區域，而`geoplot`中的`extent`引數嚴格要求經度**必須在-180到180度之間，緯度在-90到90度之間**。因此在原圖的基礎上我們進行微調，將標題移動到居中位置，具體程式碼如下：
```Python
from palettable.cartocolors.sequential import Teal_7_r
import matplotlib.font_manager as fm
from shapely.geometry import box
# 讀入世界主要河流線資料
world_river = gpd.read_file('geometry/world_rivers_dSe.geojson')
# 讀入世界海洋麵資料
world_ocean = gpd.read_file('geometry/world_ocean.shp')
# 圖層1：世界範圍背景色，基於shapely.geometry中的bbox來生成矩形向量
ax = gplt.polyplot(df=gpd.GeoDataFrame({'geometry': [box(-180, -90, 180, 90)]}),
facecolor='#000026',
edgecolor='#000026')
# 圖層2：世界海洋麵圖層
ax = gplt.polyplot(world_ocean,
facecolor='#00003a',
edgecolor='#00003a',
ax=ax)
# 圖層3：世界主要河流線圖層
ax = gplt.sankey(world_river,
scale='StrokeWeig',
hue='StrokeWeig',
scheme=mc.Quantiles(world_river['StrokeWeig'], 7),
cmap=Teal_7_r.mpl_colormap,
limits=(0.05, 0.4),
figsize=(8, 8),
extent=(-180, -90, 180, 90),
ax=ax)
# 新增標題
ax.text(0, 0, 'The Rivers of the World',
fontproperties=fm.FontProperties(fname="AlexBrush-Regular.ttf"), # 傳入Alex Brush手寫字型檔案
fontsize=28,
color=Teal_7_r.mpl_colors[-1],
horizontalalignment='center',
verticalalignment='center')
# 新增作者資訊及資料來源
ax.text(0, -15, 'Visualization by CNFeffery - Data by Natural Earth',
fontproperties=fm.FontProperties(fname="AlexBrush-Regular.ttf"),
fontsize=8,
color='#599bae',
horizontalalignment='center',
verticalalignment='center')
plt.savefig('圖13.png', dpi=600, pad_inches=0, bbox_inches='tight')
```
圖13
以上就是本文的全部內容，我將在下一篇文章中繼續與大家一起探討學習`geoplot`中更高階的繪圖API。如有疑問和意見，歡迎留言或在我的`Github`倉庫中發起*issues*與我