桑基图还能在地图上绘制？！这个可视化技巧超常用~~

本节提要：简单介绍使用geoplot来绘制地理桑基图（sankey）

---

---

前不久群里有个同学问能不能画一张漂亮的桑基图，原图找不到了，大概像下面这张。

我回答目前常用的库包不能直接绘制这样的桑基图，我错了，应该回答是目前常用的库包不能绘制这样漂亮些的桑基图。

其实geoplot库包已经内置了sankey这个命令，除了比较丑。

按照官网的说明，geoplot库包是基于matplotlib、cartopy的高级地理封装包，类似seaborn是matplotlib的高级封装包。但是真上手用起来会发现，他其实借用了很多geopandas的东西，绘图数据也以GeoDataFrame格式为主。

另外，这个库包的桑基图命令不能修改线条的宽度，所以只能通过颜色来映射数据的流向。（这就很鸡肋了）其本质是生成的带颜色映射的Line2D。其实如果不能修改线宽，还不如直接用matplotlib和cartopy硬画。

首先在虚拟环境中下载安装geoplot:

conda install -c conda-forge geoplot

sankey这个命令，需要一个GeoDataFrame来进行映射，至少需要两列数据，一列用来映射，一列用来使用geometry。

安装好库包后，导入要使用的库包：

import geoplot as gpltimport geoplot.crs as gcrsimport geopandas as gpdimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport mapclassify as mcimport numpy as npimport pandas as pdimport cartopy.io.shapereader as shpreaderfrom shapely.geometry import MultiPointextent=[108.3,109.35,29.7,30.7]plt.rcParams['font.sans-serif']=['SimHei']

接下来，读取地图文件，并转换投影使经纬度显示正常（还是费弗里大佬厉害

），给出放射中心点经纬度坐标，并随机生成用于映射的数据：

a=gpd.read_file(r'E:\2020-06-09利川市行政边界50\利川市_行政边界乡镇\利川市_行政边界.shp').to_crs('EPSG:4326')lichuan_center=(108.95,30.29)data=np.random.rand(15)

接下来使用geopandas的bounds命令，获取到每个地区的中心经纬度，当然由于行政区划的不规则性，有些点不在该行政区内部：

b=a.boundsolon=(b['minx']+b['maxx'])/2olat=(b['miny']+b['maxy'])/2

然后生成我们将要用来绘图的GeoDataFrame，通过shapely的MultiPoint功能，生成的geometry都是MULTIPOINT:

multipoints=[]for x,y in zip(olon,olat): multipoints.append(MultiPoint([lichuan_center, (x, y)]))d = {'data': data, 'geometry': multipoints}gdf=gpd.GeoDataFrame(d, crs=4326)

fig=plt.figure(figsize=(2,2),dpi=500)ax=fig.add_axes([0,0,1,1])

接下来绘制桑基图和地图以及末端的散点：

gplt.sankey(ax=ax,df=gdf,hue='data',scale='data',cmap='Reds',lw=1.5,zorder=1)gplt.polyplot(a, ax=ax, facecolor='lightgray', edgecolor='k',lw=0.5)ax.scatter(olon,olat,s=2,color='k',zorder=2)ax.set(xlim=(108.3,109.35),ylim=(29.7,30.7))ax.axis('off')

大概就是这样画出来，但我感觉实在太丑。虽然我们可以通过修改一定数量的美化参数来订正一定的形式：

但是geoplot的sankey命令最终是基于matplotlib的line2d类，这个类的线宽参数linewidth只能接受标量而不能接受可迭代的量，所以宽度是不能随每根线而变化。为了实现这种变化，我们只能定义一个函数，来绘制线宽随线值变化的桑基图，这里简单做一个事例：

import matplotlib.pyplot as pltimport numpy as npimport cartopy.crs as ccrsimport cartopy.io.shapereader as sprimport cartopy.mpl.ticker as cmtimport matplotlib.ticker as mtickerimport geopandas as gpdplt.rcParams['font.sans-serif']=['KaiTi']shp_path=r'E:\家园\利川市地图\利川.shp'fig=plt.figure(figsize=(2,2),dpi=500)ax=fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())ax.add_geometries(spr.Reader(shp_path).geometries(),crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',lw=0.5)ax.set_extent([108.3,109.35,29.7,30.7],crs=ccrs.PlateCarree())ax.set_xticks(np.arange(108.3,109.35,0.21))ax.set_yticks(np.arange(29.7,30.7,0.2))ax.xaxis.set_major_formatter(cmt.LongitudeFormatter())ax.yaxis.set_major_formatter(cmt.LatitudeFormatter())ax.tick_params(direction='in',labelsize=3,top=True,right=True,length=2,width=0.5)gl=ax.gridlines(draw_labels=False,linewidth=0.5,linestyle='--')gl.xlocator=mticker.FixedLocator(np.arange(108.3,109.35,0.21))gl.ylocator=mticker.FixedLocator(np.arange(29.7,30.7,0.2))center=(108.95,30.29)a=gpd.read_file(r'E:\2020-06-09利川市行政边界50\利川市_行政边界乡镇\利川市_行政边界.shp').to_crs('EPSG:4326')b=a.boundslons=(b['minx']+b['maxx'])/2lats=(b['miny']+b['maxy'])/2data=np.random.rand(15)def draw_sankey(ax,center,lon,lat,data):    for lo,la,d,in zip(lon,lat,data):        ax.annotate('',                    xy=(lo,la),                    xytext=center,                    size=20,                     va='center',                    ha='center',                    arrowprops=dict(connectionstyle='arc3,rad=-0.2',                                    edgecolor='none',                                    width=d/data.max()*3,headwidth=d/data.max()*5) )        ax.scatter(lo,la,c='k',s=5)        ax.text(lo,la+0.02,'{:.1f}'.format(d),fontsize=3)    return axdraw_sankey(ax,center,lons,lats,data)plt.show()

封装好的地理桑基图的绘制可定制化效果比较差，matplotlib自带的桑基命令不能和cartopy一起用。只能迂回到注释语句annotate或者arrow来画比较像的地理桑基图。不知道费弗里大佬将来会不会推出这类地图的完全geopandas的绘制方法。