如何用Python优雅地绘制中国的地图

相信很多气象专业的小伙伴在做科研的时候大概率都会遇到一个问题：如何找到中国地图的边界数据？如何绘制合规的中国地图？如何利用地图边界做裁减？如果你正在面临或者即将面临这些问题，请阅读本文，它将会极大地提高你在地图绘制上的生产力。

地图这个坑挺折磨人的，一方面我们需要使用合规的地图数据，而寻找和验证合规数据本身就是一件耗费时间和精力的事情，另一方面，拿到数据以后，我们还需要花时间去理解所拿到的数据，以shapefile格式文件为例，它是二进制格式，是人类无法直接阅读的，你需要借助一些工具来对它进行可视化，比如QGIS或者ArcGIS，而当你需要画图的时候，还需要理解shapefile这种存储格式以及相关的模块包的调用方法，最终才能正确地拿到你想要的lat/lon数据，这一整套搞下来那叫一个酸爽。

当然，我们可以在一些论坛和博客里得到很多帮助，有很多人愿意分享自己的经验，有的小伙伴会分享自己的合规的shapefile文件，有的小伙伴会分享自己利用shapefile进行剪切（白化）、画图的代码，之前我看了很多类似的文章，比如之前我还是个菜狗的时候，气象家园博主平流层的萝卜写的基于basemap的地图裁减的文章《[经验总结]Python完美白化》让我收获很大，在此表示感谢。

后来我发现，虽然写分享的人很多，但是这种分享的方式有几个的问题：1. 代码与数据分离，需要额外准备数据，如果博主不提供对应的数据那么代码一点用处也没有。2. 使用不方便，这类文章所提供的代码片段，可能需要你自己去封装，如果要在不同的项目里使用，就需要不停地copy代码。

而在我看来，解决地图绘制痛点最优雅的方式应该是直接写一个python第三方扩展包，预置数据，pip安装，极简调用。但是我似乎并没有找到有人写这样的包，既然之前没有人写，那我就来做第一个吧。

cnmaps: 让中国地图画起来更丝滑

编写cnmaps的目标是简化一切能简化的步骤，让中国区域的地图绘制的功能变得简单、直接。

目前cnmaps仅支持python>=3.6的版本，因此我们先用conda创建一个基于python3.6的干净的虚拟环境并且激活它

$ conda create -n cnmaps -y python=3.6$ conda activate cnmaps

在安装cnmaps之前，我们还需要预先安装cartopy，且建议cartopy的版本在 0.19.0 及以上，由于cartopy的依赖比较复杂，用pip安装会比较费事，所以我们用conda安装。

(cnmaps) $ conda install -c conda-forge -y cartopy==0.19.0

把上面的准备工作都做完以后，我们就可以安装cnmaps了，安装cnmaps很简单，可以直接用pip安装

(cnmaps) $ pip install cnmaps

自带地图

cnmaps想要解决的第一个痛点是：地图数据查找。

因此cnmaps是自带地图的，当你安装cnmaps以后，它就已经预先把中国国界和省界的边界数据装好了，原始数据来自高德开放API，绝对合规，可以放心食用。

让我们先来绘制一张最简单的中国地图

import cartopy.crs as ccrsimport matplotlib.pyplot as pltfrom cnmaps import get_map, draw_map
fig = plt.figure(figsize=(10,10))ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_map(get_map('中国'), color='k')draw_map(get_map('南海'), color='k')
plt.savefig('china.png', bbox_inches='tight')

我们再来绘制一张四省地图

import cartopy.crs as ccrsimport matplotlib.pyplot as pltfrom cnmaps import get_map, draw_map
PROVINCES = ['河北', '河南', '山东', '山西']
fig = plt.figure(figsize=(10, 10))fig.tight_layout()
for i, prov in enumerate(PROVINCES):    ax = fig.add_subplot(2, 2, i + 1, projection=ccrs.PlateCarree())    draw_map(get_map(prov), color='k')    plt.savefig('./prov4.png', bbox_inches='tight')

这里讲解一下，绘制地图的核心函数有两个：一个是get_map，另一个是draw_map，其中get_map传入的第一个参数是地图的名称，例如get_map('河南')，它返回的是一个MapPolygon对象，draw_map所需要接收的第一个参数就是这个MapPolygon对象，然后绘图，另外draw_map函数最终执行绘图的是plt.plot函数，dew_map除了第一个参数是地图对象以外，其他关键字参数都是传给plt.plot，注意避免传入plt.plot不认识的参数。

大家可以自行尝试感兴趣的省份来绘图，直接向get_map中传入该省的中文简写、全称、车牌的前缀以及英文名，例如get_map('山西')、get_map('山西省')、get_map('晋')和get_map('Shanxi')是等同的，都可以取到山西省的边界数据。

注意，目前cnmaps的数据尚未支持到市、县界，仅能绘制出省（直辖市、特区）级别的行政区域。

自由组合你需要的地图区域

cnmaps想要解决的第二个痛点是：地图之间的组合。

设想一下，你现在正在研究一个关于京津冀的污染的项目，需要京津冀的边界文件，这时候你要怎么做呢？可能需要先找到北京、天津和河北的shapefile文件，然后打开QGIS或者ArcGIS，去网上查找相关的教程进行操作，把三个边界合成一个，然后再输出成新的shapefile再去画图是吗？

如果用cnmaps，这一过程可以很轻松地完成。

import cartopy.crs as ccrsimport matplotlib.pyplot as pltfrom cnmaps import get_map, draw_map
jingjinji = get_map('北京') + get_map('天津') + get_map('河北')fig = plt.figure(figsize=(10, 10))ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_map(jingjinji, color='k')
plt.savefig('./jingjinji.png', bbox_inches='tight')

看到了吗？我们直接把北京、天津、河北的地图对象相加，就合并出了京津冀地区的边界，是不是很简单很直接？类似的你们可以尝试一下东三省、江浙沪、川渝...

等值线图裁减

cnmaps想要解决的第三个痛点是：繁琐的地图裁减过程。

就这个地图裁减，会在很多分享的博客文章里出现（当然包括我之前的某些文章），作者通常会贴出大量的代码细节，甚至已经细致到了对画轴笔线点的连接操作，显然这种形式会劝退很多不那么自信的读者。而某些自信的读者在粘贴了这些代码以后，又可能会因为代码对运行环境的各种依赖问题而陷入漫长的调试。

所以我认为应该有一个足够抽象、健壮、具有通用性和强内聚性的函数，让调用者完全无需知道内部细节就可以丝滑使用。因此我写了一个clip_countours_by_map函数，下面我们继续以京津冀为例子，来看一下clip_countours_by_map的效果。

import cartopy.crs as ccrsimport matplotlib.pyplot as pltimport numpy as npfrom cnmaps import get_map, draw_map, clip_contours_by_mapfrom cnmaps.sample import load_dem
lons, lats, dem = load_dem()jingjinji = get_map('北京') + get_map('天津') + get_map('河北')
fig = plt.figure(figsize=(10, 10))ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_map(jingjinji, color='k')cs = ax.contourf(lons,                 lats,                 dem,                 cmap=plt.cm.terrain,                 levels=np.linspace(-2112, dem.max(), 10))clip_contours_by_map(cs, jingjinji)ax.set_extent(jingjinji.get_extent(buffer=1))
plt.savefig('./jingjinji-cliped-dem.png', bbox_inches='tight')

clip_countours_by_map的作用是根据地图裁减等值线，它需要传入的第一个参数，是一个matplotlib.contour.QuadContourSet对象，也就是调用ax.contourf所返回的对象，第二个参数需要传入地图边界对象，也就是通过get_map以及它们之间的运算而得到的对象。

我们可以看到上述代码中又出现了一个新的方法函数get_extent，这个函数是绑定在地图对象上的，对于初学者来说，你只要知道所有get_map返回的以及它们之间运算得到的对象，都可以调用这个方法函数，你可以把它近似理解为“缩放到地图位置”。get_extent只接受一个参数buffer，不同buffer的效果大概是如下这样：

通常情况下，如果一个多边形有洞，若不加处理，在裁减时会将洞忽略掉，或者无法正确裁切，而clip_countours_by_map函数解决了类似的问题，可以正确地进行有洞多边形的裁切，例如我们把京津冀区域去掉北京（比如你在研究环京污染之类的），从而形成一个带洞的多边形再绘图。

import cartopy.crs as ccrsimport matplotlib.pyplot as pltimport numpy as npfrom cnmaps import get_map, draw_map, clip_contours_by_mapfrom cnmaps.sample import load_dem



    
lons, lats, dem = load_dem()jinji = get_map('天津') + get_map('河北')
fig = plt.figure(figsize=(10, 10))fig.tight_layout()ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_map(jinji, color='k')cs = ax.contourf(lons,                 lats,                 dem,                 cmap=plt.cm.terrain,                 levels=np.linspace(-2112, dem.max(), 10))clip_contours_by_map(cs, jinji)ax.set_extent(jinji.get_extent(buffer=1))
plt.savefig('./jinji-cliped-dem.png', bbox_inches='tight')

喏，带洞多边形的裁切就完成了。

这里我想多说两句，其实对于地图裁减，有很多种解决方案，但归根结底本质上只分为两种，一种是数据掩膜，一种是图形裁切。数据掩膜的本质是从输入端做处理，将区域外的数据做遮罩，使用这类方法的实现比较典型的就是salem，但是这类方法有一个问题，就是会有锯齿，当你数据的分辨率比较粗的时候，锯齿会非常的明显，仿佛打了码。而图形裁切方案，是从matplotlib图形对象入手，在输出端进行裁切，它的过程是先把完整的数据图整个画出来（当然也会包含插值），然后按照轮廓把所需的部分像剪纸一样“剪”出来，这样出来的图形完全不会有锯齿。

当然了，除了上述的裁减方式以外，我们还可以配合contour和clabel来绘制等值线及标签。

import cartopy.crs as ccrsimport matplotlib.pyplot as pltimport numpy as npfrom cnmaps import (get_map, draw_map, clip_contours_by_map,                    clip_clabels_by_map)from cnmaps.sample import load_dem
lons, lats, dem = load_dem('京津冀')
# 构建京津冀地图jingjinji = get_map('北京') + get_map('天津') + get_map('河北')
fig = plt.figure(figsize=(10, 10))fig.tight_layout()
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())draw_map(jingjinji, color='k')
csf = ax.contourf(lons,                  lats,                  dem,                  cmap=plt.cm.terrain,                  levels=np.linspace(-2112, 8700, 20))# 切填色等值线clip_contours_by_map(csf, jingjinji)
cs = ax.contour(lons,                lats,                dem,                colors='b',                levels=np.linspace(-2112, 8700, 20))# 切等值线clip_contours_by_map(cs, jingjinji)
clabels = ax.clabel(cs, inline=True, fmt='%i')# 用地图对象切出标签clip_clabels_by_map(clabels, jingjinji)
ax.set_extent(jingjinji.get_extent(buffer=1))
plt.savefig('./jingjinji-cliped-dem-with-clabels.png', bbox_inches='tight')

从这个例子同时也可以看出，clip_countours_by_map不仅适用于contourf，也适用于contour。

但是需要注意的一点是，由于cartopy本身在重写matplotlib的Axes类方法的时候，不知道是因为设计缺陷还是什么原因，在0.18.0版本中对于clabel方法的返回值与原版不一致，它什么也不返回，导致我们无法编辑Text对象列表，进而无法实现对clabel的裁切，好在这一问题在0.19.0及以上版本得到了修复，所以强烈建议使用cartopy>=0.19.0版本。

多种投影支持

cnmaps想要解决的第四个痛点是：投影转换。

地图投影是气象学避不开的话题，如果你研究的是行星尺度、或者高纬地区，那更是可能把各种五花八门的投影给祭出来了，所以投影问题当然是cnmaps必须要考虑的。

其实平心而论，我个人认为basemap的坐标投影转换写得很良心，使用起来简直就是无脑操作，而cartopy在坐标投影转换上的实现语法就显得非常啰嗦。但是没关系，cnmaps会帮你隐藏这些无聊的语法。

我们先来看一个例子

import cartopy.crs as ccrsimport matplotlib.pyplot as pltfrom cnmaps import get_map, draw_map, clip_contours_by_mapfrom cnmaps.sample import load_dem
lons, lats, dem = load_dem()
fig = plt.figure(figsize=(16, 12))fig.tight_layout()
china = get_map('中国')ax = fig.add_subplot(111, projection=ccrs.Robinson(central_longitude=100))cs = ax.contourf(lons,                 lats,                 dem,                 cmap=plt.cm.terrain,                 transform=ccrs.PlateCarree())clip_contours_by_map(cs, china)
draw_map(china, color='r')draw_map(get_map('南海'), color='r')ax.set_global()ax.coastlines()plt.title('Robinson')
plt.savefig('./china-clip-robinson.png', bbox_inches='tight')

对于cartopy来说，切换投影需要改的地方有两处，一是在添加子图的时候，需要把projection参数改成你想要展示的投影的cartopy定义的CRS对象，然后在调用contourf或者contour绘图的时候，把数据源所使用的原始投影的CRS对象传入transform参数即可。

其实上面描述的过程，就是cartopy给我们指定的投影转换方法（文档连接），而对于cnmaps所要修改的地方有多少呢？答案是0，也就是说只要你遵循cartopy给你指定的投影转换方法进行投影转换，那么cnmaps会自动帮你把剩下的事情通通搞定。

最后让我们来欣赏一下cnmaps在其他投影上的表现。

项目链接：https://github.com/Clarmy/cnmaps

欢迎提交Issue和PR

也欢迎访问我重新开张的个人博客clarmy.net😈