点击蓝字关注我，有干货领取！

本文对应代码和数据已上传至我的Github仓库：
https://github.com/CNFeffery/DataScienceStudyNotes^[1]

已发布：

Python地信专题 | 基于geopandas的空间数据分析—数据结构篇

Python地信专题 | 基于geopandas的空间数据分析-坐标参考系篇

Python地信专题 | 基于geopandas的空间数据分析-文件IO篇

Python地信专题 | 基于geopandas玩转地图可视化

1 简介

通过前面的文章，我们已经对geopandas中的数据结构、坐标参考系、文件IO以及基础可视化有了较为深入的学习。

其中在基础可视化那篇文章中我们提到了分层设色地图，可以对与多边形关联的数值属性进行分层，并分别映射不同的填充颜色。

但只是开了个头举了个简单的例子，实际数据可视化过程中的分层设色有一套策略方法。

作为基于geopandas的空间数据分析系列文章的第五篇，通过本文你将会学习到基于geopandas和机器学习的分层设色。

2 基于geopandas的分层设色

地区分布图（Choropleth maps，又叫面量图）作为可能是最常见的一种地理可视化方法。

其核心是对某个与矢量面关联的数值序列进行有意义的分层，并为这些分层选择合适美观的色彩，最后完成对地图的着色。

优点是美观且直观，即使对地理信息一窍不通的人，也能通过颜色区分出不同面之间的同质性与异质性：

但同样地，如果对数据分层采取的方法有失严谨没有很好的遵循数据特点，会很容易让看到图的人产生出不正确的判断。

下面我们按照先分层，后设色的顺序进行介绍。

2.1 基于mapclassify的数据分层

上一篇文章中我们提到过,，在geopandas.GeoDataFrame.plot()中，参数scheme对应的数据分层是基于第三方库mapclassify实现的。

因此要想对geopandas中的数据分层有深入的了解，我们就得先来了解一下mapclassify中的各种数据分层算法。

用到的数据是系列文章前几期使用地滚瓜烂熟的新冠肺炎疫情数据，数据处理过程同上一篇文章，这里不再解释：

2.1.1 BoxPlot

BoxPlot即箱线图，是统计学中使用到的一种方法：

对个数为  观测数据从小到大进行排序，分别得到位置处于  、 以及  的观测值，称为  、 以及  （即第一四位数、中位数和第三四分位数）。

并定义  为  ，以  为下限，以  为上限，将小于下限或大于上限的观测值作为离群异常值。

最后用图像的形式表达上述计算结果，如图2的上图，而图2的下图对应着概率估计。

可以看出，箱线图法实际上是基于概率估计的一种异常值剔除方法，因为离群值只有  的概率会出现，即如果你想要找出数据中的异常高低值，BoxPlot是不错的选择：

在mapclassify中我们使用BoxPlot()来为数据实现箱线图分层：

import mapclassify as mc

# 对各省2020-03-08对应的累计确诊数量进行分层
bp = mc.BoxPlot(temp['province_confirmedCount'])
# 查看数据分层结果
bp

可以看出通过箱线图法将数据分成了五类，其中异常值只有1个即为湖北省。

下面我们配合geopandas来对上述结果进行可视化，和上一篇文章一样，按照省级单位名称连接我们的疫情数据与矢量数据：

接着对其进行可视化，在上一篇文章图28的基础上，将scheme参数改为BoxPlot，又因为箱线图可以看作无监督问题，故分层数量k在这里无效，删去：

fig, ax = plt.subplots(figsize=(10, 10))

ax = data_with_geometry.to_crs(albers_proj).plot(ax=ax,
                                                 column='province_confirmedCount',
                                                 cmap='Reds',
                                                 missing_kwds={
                                                     "color": "lightgrey",
                                                     "edgecolor": "black",
                                                     "hatch": "",
                                                     "label": "缺失值"
                                                 },
                                                 legend=True,
                                                 scheme='BoxPlot',
                                                 legend_kwds={
                                                     'loc': 'lower left',
                                                     'title': '确诊数量分级',
                                                     'shadow': True
                                                 })

ax = nine_lines.geometry.to_crs(albers_proj).plot(ax=ax,
                                                  edgecolor='grey',
                                                  linewidth=3,
                                                  alpha=0.4)

ax.axis('off')
plt.suptitle('新型冠状肺炎累计确诊数量地区分布', fontsize=24) # 添加最高级别标题
plt.tight_layout(pad=4.5) # 调整不同标题之间间距
ax.text(-2800000, 1300000, '* 原始数据来源：丁香园，n其中台湾及香港数据缺失') # 添加数据说明

fig.savefig('图6.png', dpi=300)

咋看起来没问题，但是如果你仔细观察左下角的图例会发现前两行范围颜色是重复的，且数值范围是错乱的。

这是geopandas.GeoDataFrame.plot()中涉及箱线图法的一个小bug，遇到这种问题不用慌。

如果你在上一篇文章中去我的Github仓库查看过创作图29对应的代码，一定会想到既然geopandas自身有bug，那我们用matplotlib中的mpatches和legend自定义图例就可以啦。

而为了自定义的图例色彩与geopandas映射出的保持一致，我们需要额外使用到matplotlib中的get_cmap(cmap)来制作可独立导出颜色的cmap方案实例。

譬如我们这里是Reds，就需要按照前面bp的有记录数量的分层结果，从Reds中产生同样5个档次的颜色，具体操作过程如下：

import matplotlib.patches as mpatches

fig, ax = plt.subplots(figsize=(10, 10))

ax = data_with_geometry.to_crs(albers_proj).plot(ax=ax,
                                                 column='province_confirmedCount',
                                                 cmap='Reds',
                                                 missing_kwds={
                                                     "color": "lightgrey",
                                                     "edgecolor": "black",
                                                     "hatch": "",
                                                     "label": "缺失值"
                                                 },
                                                 scheme='BoxPlot')

handles, labels = ax.get_legend_handles_labels() #get existing legend item handles and labels

ax = nine_lines.geometry.to_crs(albers_proj).plot(ax=ax,
                                                  edgecolor='grey',
                                                  linewidth=3,
                                                  alpha=0.4)

# 实例化cmap方案
cmap = plt.get_cmap('Reds')

# 得到mapclassify中BoxPlot的数据分层点
bp = mc.BoxPlot(temp['province_confirmedCount'])
bins = bp.bins

# 制作图例映射对象列表
LegendElement = [mpatches.Patch(facecolor=cmap(_*0.25), label=f'{int(max(bins[_], 0))} - {int(bins[_+1])}')
                 for _ in range(5)] + 
                [mpatches.Patch(facecolor='lightgrey', edgecolor='black', hatch='', label='缺失值')]

# 将制作好的图例映射对象列表导入legend()中，并配置相关参数
ax.legend(handles = LegendElement, loc='lower left', fontsize=10, title='确诊数量分级', shadow=True, borderpad=0.6)

ax.axis('off')
plt.suptitle('新型冠状肺炎累计确诊数量地区分布(截至2020年03月04日)', fontsize=24) # 添加最高级别标题
plt.title('数据分层方法：BoxPlot', fontsize=18)
plt.tight_layout(pad=4.5) # 调整不同标题之间间距
ax.text(-2900000, 1250000, '* 原始数据来源：丁香园，n其中台湾及香港数据缺失') # 添加数据说明

fig.savefig('图7.png', dpi=300)

可以看到，通过自定义图例的方式，虽然麻烦了一点，但是我们不仅修复了图例的bug，还为其添加了更加完善的细节。

如图形修改为矩形，范围修改为整数。

2.1.2 EqualInterval

EqualInterval即等间距，是最简单的一种分层方法。

它在原数据最小值与最大值间以等间距的方式划分出k个层次，mapclassify中对应等间距法的类为EqualInterval()：

bp = mc.EqualInterval(temp['province_confirmedCount'])
# 查看数据分层结果
bp

可以看到对于分布非常不均匀的新冠肺炎确诊数量数据来说，这种方法表现得十分糟糕，中间三个类都没有记录落入。

如果使用这种方法强行绘图，效果就会类似上一篇文章中地区分布图部分。

最开始那个糟糕的效果那样只有湖北一个地方是最深的暗红色，而其他地方皆为最淡的色阶，这里就不重复演示。

2.1.3 FisherJenks

在了解mapclassify中的FisherJenks之前，我们先来了解一下什么是Jenks Natural Breaks：

• Jenks Natural Breaks

Jenks Natural Breaks旨在为1维数据计算合适的划分点，使得不同组之间的差距尽可能大的同时组内差距尽可能小。

其思路非常简单，举一个简单的例子进行说明：

对于一组待分割的序列  ，现在需要为其找到将原始数据分为  部分的方法。

那么实际上就有  、 以及  这三种切分方法，现定义sum of squared deviations for array mean（简称SDAM）：