前言

最早笔者用 Python 画中国地图时，会准备 bou2_4p.shp 文件，然后封装一个读取 shapefile 并添加到 GeoAxes 上的函数，别的项目要用时就把数据和函数复制粘贴过去。Cartopy 系列：从入门到放弃 里就是这么做的。

后来工作中用到了 Clarmy 开发的 cnmaps 包，只用两行就能快速绘制地图，非常方便。同时萌生了自己实现一个功能类似的包的想法，遂开发出了 frykit。

前言

裁剪或者说白化，就是让填色图只显示在多边形里面，不显示在多边形外面，例如只显示 GeoAxes.contourf 在中国境内的结果。实现方法为：

from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
from cartopy.io.shapereader import Reader

reader = Reader(filepath)
geom = next(reader.geometries())
reader.close()

cf = ax.contourf(X, Y, Z, transform=crs)
geom = ax.projection.project_geometry(geom, crs)
path = Path.make_compound_path(*geos_to_path(geom))
for col in cf.collections:
    col.set_clip_path(path, ax.transData)

• 将 crs 坐标系上的多边形对象变换到 data 坐标系上。
• 利用 geos_to_path 和 make_compound_path 将变换后的多边形转为 Path 对象。
• 对 QuadContourSet.collections 里的每个成员调用 set_clip_path 方法，并且指定 data 坐标系。

二维平面上一系列点的坐标由 x 和 y 数组描述，同时准备一个形状相同的 mask 数组。若第 i 个点落入了平面上一个多边形的内部，则令 mask[i] = True；若在多边形外，则令 mask[i] = False。由此得到的 mask 数组即掩膜（mask）数组，它可以作为布尔索引分出多边形内外的点

x_in, y_in = x[mask], y[mask]
x_out, y_out = x[mask], y[mask]

它可以作为掩膜，掩盖多边形范围外的值——即把外面的值设为 NaN，以便进行后续的计算

z[~mask] = np.nan
z_mean = np.nanmean(z)

下图展示了两个应用：左小图的多边形是一个中心带洞的正方形，给定一系列散点的坐标，计算出掩膜后可以把多边形内的散点画成红色，多边形外的散点画成蓝色；右小图的多边形是中国全域，给定 (50, 50) 形状的经纬度网格，计算出掩膜后用橙色画出掩膜为 True 的部分，这张掩膜之后可以用来处理网格上的其它变量。

本文的目的是介绍如何用 Python 制作掩膜数组，并尽量优化其运行时间。从 shapefile 中读取中国国界并转化为 Shapely 中的多边形对象等操作，已经在博文 Cartopy 系列：探索 shapefile 中详细介绍过了，本文是对其的一个补充。

前言

Cartopy 可以通过 feature 模块向地图添加国界 BORDER 和省界 STATES，因其底层采用的 Natural Earth 地图数据并不符合我国的政治主张，所以我们经常需要自备 shapefile 文件来画中国省界，以下面的代码为例

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

extents = [70, 140, 0, 60]
crs = ccrs.PlateCarree()
fig = plt.figure()
ax = fig.add_subplot(111, projection=crs)
ax.set_extent(extents, crs)

filepath = './data/bou2_4/bou2_4p.shp'
reader = shpreader.Reader(filepath)
geoms = reader.geometries()
ax.add_geometries(geoms, crs, lw=0.5, fc='none')
reader.close()

plt.show()

图就不放了，这段代码足以应付大部分需要画省界的情况。然而我在无脑粘贴代码的过程中逐渐产生了疑惑：为什么 shapefile 会由三个文件组成？省界是以何种形式存储在文件中？Cartopy 和 Matplotlib 又是怎样将省界画出来的？调查一番源码后总结出了这段代码底层实现的流程：

• 利用 PyShp 包读取 shapefile 文件中的每个形状。
• 利用 Shapely 包将形状转换为几何对象。
• 利用 Cartopy 包将几何对象投影到地图所在的坐标系上。
• 用投影后的坐标构造 Matplotlib 的 Path 对象，最后画在地图上。

本文的目的即是从头到尾解说一下这段流程，希望加深对 shapefile 格式，Matplotlib 和 Cartopy 包的理解。令人意外的是，随着探索的不断深入，我发现自己自然而然地学会了如何实现省份填色、省份合并，地图白化等，以前看起来十分困难的操作。本文也会一并介绍这些应用。

前言

Matplotlib 中用箭头表示风场或电磁场等矢量场时需要用到 quiver 方法，据字典，quiver 一词的意思是颤动、颤抖或箭袋，貌似也就最后一个意思跟箭头搭得上边。相比于其它画图方法，quiver 的参数又多又容易混淆，所以本文将以图解的方式逐一介绍。这些参数按功能可分为三种：控制箭头位置和数值的、控制箭头长度和角度的，以及控制箭头尺寸和形状的。下面会按照这个分组顺序来解说。本文代码基于 Matplotlib 3.3.4。

前言

几年前曾写过 Cartopy 系列：从入门到放弃，但现在来看还是遗漏了不少细节，比如初学者可能会遇到以下问题

• 经度是用 [-180°, 180°] 还是 [0°, 360°] 范围？
• 为什么有时候设置的刻度显示不全？
• 怎么截取跨越地图边界的区域，画图又怎么跨越边界？

本文将会用坐标变换的思想来解答以上问题，希望能给读者一些实用的启示。本来应该把这些内容写到入门教程里的，但可能会太长，所以现在单独成篇。文中的讨论主要针对最常用的 Plate Carrée 投影，其它投影需要读者自己测试。代码基于 Cartopy 0.18.0，虽然现在已经更新到 0.20.0 了，但基本思想是一致的。

前言

Cartopy 中的 Plate Carrée 投影使用方便，但在展示中国地图时会使中国的形状显得很瘪，与之相比，Lambert 投影的效果会更加美观，下图显示了两种投影的差异

所以本文将会介绍如何在 Cartopy 中实现 Lambert 投影，并为地图添上合适的刻度。文中 Cartopy 的版本是 0.18.0。

简介

常用的地图可视化的编程工具有 MATLAB、IDL、GrADS、GMT、NCL 等。我之前一直使用的是脚本语言 NCL，易用性不错，画地图的效果也很好。然而 2019 年初，NCAR 宣布 NCL 将停止更新，并会在日后转为 Python 的绘图包。于是我开始考虑转投 Python，同时觉得在 Python 环境下如果还是用 PyNGL 那一套语法的话，未免有些换汤不换药。因此我选择用 Python 环境下专有的 Cartopy 包来画地图。

此前 Python 最常用的地图包是 Basemap，然而它将于 2020 年被弃用，官方推荐使用 Cartopy 包作为替代。Cartopy 是英国气象局开发的地图绘图包，实现了 Basemap 的大部分功能，还可以通过 Matplotlib 的 API 实现丰富的自定义效果。

本文将会从一个 NCL 转 Python 的入门者的角度，介绍如何安装 Cartopy，如何绘制地图，并实现一些常用的效果。代码基于 0.18.0 版本的 Cartopy。