如何更优雅地绘制/白化中国的地图

clarmy

更新：cnmaps的版本已经更新，下文中的get_map函数在新版本中已经弃用，新版本的使用指南：https://cnmaps-doc.readthedocs.io/zh_CN/latest/index.html

                               
登录/注册后可看大图

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


地图这个坑挺折磨人的，一方面我们需要使用合规的地图数据，而寻找和验证合规数据本身就是一件耗费时间和精力的事情，另一方面，拿到数据以后，我们还需要花时间去理解所拿到的数据，以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==0.2.1



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

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

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

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. from cnmaps import get_map, draw_map
   
4. 
   
5. fig = plt.figure(figsize=(10,10))
   
6. ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
   
7. draw_map(get_map('中国'), color='k')
   
8. draw_map(get_map('南海'), color='k')
   
9. plt.savefig('china.png', bbox_inches='tight')

复制代码

自由组合你需要的地图区域cnmaps想要解决的第二个痛点是：地图之间的组合。

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

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

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. from cnmaps import get_map, draw_map
   
4. 
   
5. jingjinji = get_map('北京') + get_map('天津') + get_map('河北')
   
6. 
   
7. fig = plt.figure(figsize=(10, 10))
   
8. 
   
9. ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
   
10. draw_map(jingjinji, color='k')
    
11. 
    
12. plt.savefig('./jingjinji.png', bbox_inches='tight')

复制代码

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

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

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






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

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. import numpy as np
   
4. from cnmaps import get_map, draw_map, clip_contours_by_map
   
5. from cnmaps.sample import load_dem
   
6. 
   
7. lons, lats, dem = load_dem()
   
8. 
   
9. jingjinji = get_map('北京') + get_map('天津') + get_map('河北')
   
10. 
    
11. fig = plt.figure(figsize=(10, 10))
    
12. 
    
13. ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    
14. draw_map(jingjinji, color='k')
    
15. 
    
16. cs = ax.contourf(lons,
    
17.                  lats,
    
18.                  dem,
    
19.                  cmap=plt.cm.terrain,
    
20.                  levels=np.linspace(-2112, dem.max(), 10))
    
21. clip_contours_by_map(cs, jingjinji)
    
22. ax.set_extent(jingjinji.get_extent(buffer=1))
    
23. 
    
24. plt.savefig('./jingjinji-cliped-dem.png', bbox_inches='tight')

复制代码

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

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. import numpy as np
   
4. from cnmaps import get_map, draw_map, clip_contours_by_map
   
5. from cnmaps.sample import load_dem
   
6. 
   
7. lons, lats, dem = load_dem()
   
8. 
   
9. jinji = get_map('天津') + get_map('河北')
   
10. 
    
11. fig = plt.figure(figsize=(10, 10))
    
12. fig.tight_layout()
    
13. 
    
14. ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    
15. draw_map(jinji, color='k')
    
16. 
    
17. cs = ax.contourf(lons,
    
18.                  lats,
    
19.                  dem,
    
20.                  cmap=plt.cm.terrain,
    
21.                  levels=np.linspace(-2112, dem.max(), 10))
    
22. clip_contours_by_map(cs, jinji)
    
23. ax.set_extent(jinji.get_extent(buffer=1))
    
24. 
    
25. plt.savefig('./jinji-cliped-dem.png', bbox_inches='tight')

复制代码

喏，带洞多边形的裁切就完成了。
当然了，除了上述的裁减方式以外，我们还可以配合contour
和clabel
来绘制等值线及标签。

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. import numpy as np
   
4. from cnmaps import (get_map, draw_map, clip_contours_by_map,
   
5.                     clip_clabels_by_map)
   
6. from cnmaps.sample import load_dem
   
7. 
   
8. lons, lats, dem = load_dem('京津冀')
   
9. 
   
10. # 构建京津冀地图
    
11. jingjinji = get_map('北京') + get_map('天津') + get_map('河北')
    
12. 
    
13. fig = plt.figure(figsize=(10, 10))
    
14. fig.tight_layout()
    
15. 
    
16. ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
    
17. draw_map(jingjinji, color='k')
    
18. 
    
19. csf = ax.contourf(lons,
    
20.                   lats,
    
21.                   dem,
    
22.                   cmap=plt.cm.terrain,
    
23.                   levels=np.linspace(-2112, 8700, 20))
    
24. # 切填色等值线
    
25. clip_contours_by_map(csf, jingjinji)
    
26. 
    
27. cs = ax.contour(lons,
    
28.                 lats,
    
29.                 dem,
    
30.                 colors='b',
    
31.                 levels=np.linspace(-2112, 8700, 20))
    
32. # 切等值线
    
33. clip_contours_by_map(cs, jingjinji)
    
34. 
    
35. clabels = ax.clabel(cs, inline=True, fmt='%i')
    
36. # 用地图对象切出标签
    
37. clip_clabels_by_map(clabels, jingjinji)
    
38. 
    
39. ax.set_extent(jingjinji.get_extent(buffer=1))
    
40. 
    
41. plt.savefig('./jingjinji-cliped-dem-with-clabels.png', bbox_inches='tight')

复制代码

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

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

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

我们先来看一个例子

1. import cartopy.crs as ccrs
   
2. import matplotlib.pyplot as plt
   
3. from cnmaps import get_map, draw_map, clip_contours_by_map
   
4. from cnmaps.sample import load_dem
   
5. 
   
6. lons, lats, dem = load_dem()
   
7. 
   
8. fig = plt.figure(figsize=(16, 12))
   
9. fig.tight_layout()
   
10. 
    
11. china = get_map('中国')
    
12. ax = fig.add_subplot(111, projection=ccrs.Robinson(central_longitude=100))
    
13. cs = ax.contourf(lons,
    
14.                  lats,
    
15.                  dem,
    
16.                  cmap=plt.cm.terrain,
    
17.                  transform=ccrs.PlateCarree())
    
18. clip_contours_by_map(cs, china)
    
19. 
    
20. draw_map(china, color='r')
    
21. draw_map(get_map('南海'), color='r')
    
22. ax.set_global()
    
23. ax.coastlines()
    
24. plt.title('Robinson')
    
25. 
    
26. plt.savefig('./china-clip-robinson.png', bbox_inches='tight')

复制代码

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

项目链接：https://github.com/Clarmy/cnmaps
欢迎提交Issue和PR
博客原文: http://www.clarmy.net/2022/02/12 ... teelly-with-python/

river

大神啊，必须支持一下！！！！

yqcxm

牛啊，膜拜，市县的情况呢

lovepipi628

感谢您的经验和学习分享~！

qinzhao

谢谢大神,膜拜

clarmy

yqcxm 发表于 2022-3-11 08:34
牛啊，膜拜，市县的情况呢

市县目前还没有支持，以后会有计划加上

txxrk

比maskout省劲，但是还不能做市县

clarmy

txxrk 发表于 2022-3-23 14:18
比maskout省劲，但是还不能做市县

现在可以做县市了，可以参考我的新博客文章：《cnmaps新升级:区县级数据现已加入豪华午餐》

myj0219

6666666666666666666666666

Rossby.

我的cnmaps一直安装失败，报错Timeouterror，请问该怎么解决