地理桑基图的绘制（基础知识，python）

最近年中总结，学习一下基础知识

我们的目标实现下面这张图

image.png

导入模块

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]

接下来，读取地图文件，并转换投影使经纬度显示正常，给出放射中心点经纬度坐标，并随机生成用于映射的数据：

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)

image.png

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')

image.png image.png

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

image.png image.png

但是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()

image.png image.png

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

个人爱好，看到了试一下，也许以后用得到

基础知识，这个基础知识大家知道就好，不用深入研究