(資料科學學習手札83)基於geopandas的空間資料分析——geoplot篇(下)

費弗裡發表於2020-05-16

本文示例程式碼、資料及檔案已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

  在上一篇文章中我們詳細學習了geoplot中較為基礎的三種繪圖API:pointplot()polyplot()以及webmap(),而本文將會承接上文的內容,對geoplot中較為實用的幾種高階繪圖API進行介紹。

(資料科學學習手札83)基於geopandas的空間資料分析——geoplot篇(下)
圖1

  本文是基於geopandas的空間資料分析系列文章的第7篇,通過本文你將學習geoplot中的高階繪圖API。

2 geoplot進階

  上一篇文章中的pointplot()polyplot以及webmap()幫助我們解決了在繪製散點、基礎面以及新增線上地圖底圖的問題,為了製作出資訊量更豐富的視覺化作品,我們需要更強的操縱向量資料與對映值的能力,geoplot為我們封裝好了幾種常見的高階視覺化API。

2.1 Choropleth

  Choropleth圖又稱作地區分佈圖面量圖,我們在系列之前的深入淺出分層設色篇中介紹過其原理及geopandas實現,可以通過將指標值對映到面資料上,以實現對指標值地區分佈的視覺化。

  在geoplot中我們可以通過choropleth()來快速繪製地區分佈圖,其主要引數如下:

df:傳入對應的GeoDataFrame物件

projection:用於指定投影座標系,傳入geoplot.crs中的物件

hue:傳入對應df中指定列名或外部序列資料,用於對映面的顏色,預設為None即不進行設色

cmap:和matplotlib中的cmap使用方式一致,用於控制色彩對映方案

alpha:控制全域性色彩透明度

scheme:作用類似geopandas中的scheme引數,用於控制分層設色,詳見本系列文章的分層設色篇,但不同的是在geoplot0.4.0版本之後此引數不再搭配分層數量k共同使用,而是更新為傳入mapclassify分段結果物件,下文中會做具體演示

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)

axmatplotlib座標軸物件,如果需要在同一個座標軸內疊加多個圖層就需要用這個引數傳入先前待疊加的ax

hatch:控制填充陰影紋路,詳情見本系列文章前作基礎視覺化篇圖7

edgecolor:控制多邊形輪廓顏色

linewidth:控制多邊形輪廓線型

  下面我們通過實際的例子來學習geoplot.choropleth的用法,這裡使用到的資料為美國新型冠狀肺炎各州病例數分佈,對應日期為2020年5月14日,來自Github倉庫:https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_covid_19_daily_reports_us;使用到的美國本土各州向量面資料contiguous-usa.geojson已上傳到文章開頭對應的Github倉庫中:

(資料科學學習手札83)基於geopandas的空間資料分析——geoplot篇(下)
圖2
(資料科學學習手札83)基於geopandas的空間資料分析——geoplot篇(下)
圖3

  首先我們將兩張表中各自對應的州名資料作為鍵進行連線(注意pd.merge返回的結果型別為DataFrame,需要轉換回GeoDataFrame):

# 按照州名列進行連線
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萬,遠遠超過其他州,所以這裡作為單獨的圖層進行陰影填充以突出其嚴重程度:

# 圖層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')
(資料科學學習手札83)基於geopandas的空間資料分析——geoplot篇(下)
圖4

  這樣我們就得到了圖4,需要注意的是,geoplot.choropleth()只能繪製地區分佈圖,傳入面資料後hue引數必須指定對應對映列,否則會報錯,因此這裡我們疊加紐約州和新澤西州單獨面圖層時使用的是polyplot()

2.2 Kdeplot

  geoplot中的kdeplot()對應核密度圖,其基於seaborn中的kdeplot(),通過對向量點資料分佈計算核密度估計,從而對點資料進行視覺化,可用來展示點資料的空間分佈情況,其主要引數如下:

df:傳入對應的存放點資料的GeoDataFrame物件

projection:用於指定投影座標系,傳入geoplot.crs中的物件

cmap:和matplotlib中的cmap使用方式一致,用於控制色彩對映方案

clipGeoSeries型,用於為初始生成的核密度影像進行蒙版裁切,下文會舉例說明

extent:元組型,用於傳入左下角和右上角經緯度資訊來設定地圖空間範圍,格式為(min_longitude, min_latitude, max_longitude, max_latitude)

figsize:元組型,用於控制畫幅大小,格式為(x, y)

axmatplotlib座標軸物件,如果需要在同一個座標軸內疊加多個圖層就需要用這個引數傳入先前待疊加的ax

shade:bool型,當設定為False時只有等值線被繪製出,當設定為True時會繪製核密度填充

shade_lowest:bool型,控制是否對概率密度最低的層次進行填充,下文會舉例說明

n_levels:int型,控制等值線數量,即按照概率密度對空間進行均勻劃分的數量

  下面我們回到上一篇文章開頭的例子——紐約車禍記錄資料,在其他引數均為預設的情況下,呼叫kdeplot對車禍記錄點資料的空間分佈進行視覺化:

# 圖層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
圖5

  可以看到,在kdeplot()主要引數均為預設值的情況下,我們得到了點資料空間分佈的概率估計結果及其等高線,譬如圖中比較明顯能看到的兩個點分佈較為密集的中心,下面我們調整n_levles引數到比較大的數字:

# 圖層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
圖6

  可以看到在增大n_levels引數後,圖中等值線的數量隨之增加,下面我們設定shade=True

# 圖層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
圖7

  這時影像等值線間得到相應顏色的填充,使得點分佈中心看起來更加明顯,再新增引數shade_lowest=True,即可對空白區域進行填充:

圖8
圖8

  隨之而來的問題是整幅影像都被填充,為了裁切出核密度影像的地區輪廓,將底層行政區面資料作為clip的引數傳入,便得到理想的效果:

圖9
圖9

2.3 Sankey

  桑基圖專門用於表現不同物件之間某個指標量的流動情況,譬如最常見的航線流向情況,其本質是對線資料進行視覺化,並將指標值對映到線的色彩或粗細水平上,而geoplot中的sankey()可以用來繪製這種圖,尷尬的是sankey()繪製出的OD流向圖實在太醜,但sankey()中將數值對映到線資料色彩和粗細的特性可以用來進行與流量相關的視覺化,其主要引數如下:

df:傳入對應的GeoDataFrame物件

projection:用於指定投影座標系,傳入geoplot.crs中的物件

hue:傳入對應df中指定列名或外部序列資料,用於對映線的顏色,預設為None即不進行設色

cmap:和matplotlib中的cmap使用方式一致,用於控制色彩對映方案

alpha:控制全域性色彩透明度

scheme:作用類似geopandas中的scheme引數,用於控制分層設色,詳見本系列文章的分層設色篇,但不同的是在geoplot0.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)

axmatplotlib座標軸物件,如果需要在同一個座標軸內疊加多個圖層就需要用這個引數傳入先前待疊加的ax

  下面我們以2015年華盛頓街道路網日平均交通流量資料為例,其中每個要素均為線要素,aadt代表日均流量:

圖10
圖10

  我們將其流量列對映到線的粗細程度和顏色上來,為了美觀起見我們選擇系列文章分層設色篇palettableSunsetDark作為配色方案:

# 選擇配色方案為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
圖11

2.4 geoplot中的座標參考系

  geoplot中的座標參考系與geopandas中管理起來的方式截然不同,因為geopandas基於pyproj管理座標參考系,而geoplot中的crs子模組來源於cartopy,這一點我跟geoplot的主要開發者聊過,他表示geoplot暫時不支援geopandas中那樣自定義任意投影或使用EPSG投影,而是內建了一系列常用的投影,譬如我們上文中繪製美國區域時頻繁使用到的AlbersEqualArea()即之前我們在geopandas中通過proj4自定義的阿爾伯斯等面積投影,其他常見投影譬如Web MercatorRobinson,或者直接繪製球體地圖,如本文開頭的圖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
圖12

  針對其河流寬度方面的視覺化,我們基於上文中的sankey()來實現,由於原圖中南極洲區域實際上是誇大了的,其R原始碼中設定的緯度範圍達到了-110度,這是原作者為了放得下標題內容,所以在影像下部區域虛構了一篇區域,而geoplot中的extent引數嚴格要求經度必須在-180到180度之間,緯度在-90到90度之間。因此在原圖的基礎上我們進行微調,將標題移動到居中位置,具體程式碼如下:

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
圖13

  以上就是本文的全部內容,我將在下一篇文章中繼續與大家一起探討學習geoplot中更高階的繪圖API。如有疑問和意見,歡迎留言或在我的Github倉庫中發起issues與我交流。

相關文章