python可视化 | 地理桑基图的绘制方法

郭好奇同学

发布于 2021-05-28 17:21:01

1.7K0

发布于 2021-05-28 17:21:01

文章被收录于专栏：好奇心Log

本节提要：简单介绍使用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 gplt
import geoplot.crs as gcrs
import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import mapclassify as mc
import numpy as np
import pandas as pd
import cartopy.io.shapereader as shpreader
from shapely.geometry import MultiPoint
extent=[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.bounds
olon=(b['minx']+b['maxx'])/2
olat=(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 plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as spr
import cartopy.mpl.ticker as cmt
import matplotlib.ticker as mticker
import geopandas as gpd
plt.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.bounds
lons=(b['minx']+b['maxx'])/2
lats=(b['miny']+b['maxy'])/2
data=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 ax
draw_sankey(ax,center,lons,lats,data)
plt.show()