前言

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

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

前言

人眼可见色域在色度图中表现为彩色的马蹄形，单色光（monochromatic light）的颜色对应于马蹄的弧形边界。本文想将单色光的颜色按波长线性增大的顺序一字排开，用类似彩虹渐变图的形式展示单色光光谱。用 Python 的 Matplotlib 包来实现的话，很快就能决定画图思路：

1. 读取 XYZ 颜色匹配函数（CMF）作为 XYZ 三刺激值。
2. XYZ 变换为 sRGB，接着做 gamma 校正。
3. 用 RGB 数组构造 ListedColormap 对象，用 plt.colorbar 画出。

RGB 要求范围在 $[0, 1]$，但 CMF 直接计算出的 RGB 既有负数分量，也有大于 1 的分量，所以必须采用一种方法处理范围外的分量。最后的画图效果会因处理方法的不同产生很大差别，例如下图的三条光谱：

就采取了不同的处理方式，因此在发色、颜色过渡，和亮度表现上都大有不同。本文将尝试实现不同的效果并加以分析。完整代码和相关数据见 我的 Github 仓库。

前言

1920 年代末 Wright 和 Guild 的颜色匹配实验发展出了用红绿蓝三基色（primaries）定量表示所有人眼可见颜色的 CIE RGB 色彩空间，1931 年国际照明委员会（CIE）通过对 CIE RGB 色彩空间做线性变换得到了 CIE XYZ 色彩空间。XYZ 空间里的人眼可见色域（gamut of human vision）是一块从原点出发，向无限远处不断延伸的立体区域。将这块区域投影到 $X + Y + Z = 1$ 的平面上，就能画出方便展示的 CIE 1931 色度图（chromaticity diagram）（图自 维基）：

前言

裁剪或者说白化，就是让填色图只显示在多边形里面，不显示在多边形外面，例如只显示 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 坐标系。

前言

CALIPSO 卫星的 L2 VFM（Vertical Feature Mask）产品根据激光的后向散射和消光信息，将激光通过的各高度层分类为云或气溶胶。该产品在现实中的表现如下图所示：卫星一边在轨道上移动一边向地面发射激光脉冲，相当于在地面上缓缓拉开一幅“画卷”，VFM 描述了“画卷”上云和气溶胶的分布和分类情况。

处理 VFM 产品的难点在于：

• VFM 数组呈 (N, 5515) 的形状，N 表示卫星移动时产生了 N 次观测，但 5515 并非表示有 5515 层高度，而是三种水平和垂直分辨率都不同的数据摊平成了长 5515 的数组。因此处理数据时需要参照文档的说明对 5515 进行变形。
• 文件中的经纬度和时间与 5515 的对应关系。时间数组需要解析成可用的格式。
• 每个 range bin 的分类结果编码到了 16 位的无符号短整型的每个比特上，需要按位解码。
• 网上现成的代码偏少。

网上能找到的代码有：

• CALIOPmatlab：以前 VFM 的在线文档里是给出过 MATLAB 和 IDL 的代码的，但现在链接消失了。这个仓库提供了民间改进后 MATLAB 代码。
• HDF-EOS COMPREHENSIVE EXAMPLES：HDF-EOS 网站的示例，简单易理解。
• MeteoInfo examples: CALIPSO data：基于 MeteoInfo 的代码，还有其它产品的例子。
• Visualization of CALIPSO (VOCAL)：CALIPSO 官方基于 Python 2 的可视化工具。
• 星载激光雷达CALIPSO-VFM产品数据读取与显示：MATLAB 代码的讲解。

笔者也曾写过两次教程：

• NCL绘制CALIPSO L2 VFM图像：写得很烂，作图部分可能存在问题。
• Python 绘制 CALIPSO L2 VFM 产品

本文是对旧教程的翻新，会对 VFM 数据的结构进行更多解释，对代码也进行了更新。本文使用 pyhdf 读取 HDF4 文件，用 Matplotlib 3.6.2 画图。为了方便画图，用了一些自制的函数（frykit）。虽然基于 Python，但希望能给使用其它语言的读者提供一点思路。

完整代码已放入仓库 calipso-vfm-visualization。

前言

Matplotlib 中画折线图用 ax.plot(x, y)，当横坐标 x 是时间数组时，例如 datetime 或 np.datetime64 构成的列表，x 和 y 的组合即一条时间序列。Matplotlib 能直接画出时间序列，并自动设置刻度。下面以一条长三年的气温时间序列为例：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv('test.csv', index_col=0, parse_dates=True)
series = df.loc['2012':'2014', 'T']

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(series.index, series)
ax.set_ylabel('Temperature (℃)')

print(ax.xaxis.get_major_locator())
print(ax.xaxis.get_major_formatter())

<matplotlib.dates.AutoDateLocator object at 0x000001AC6BF89A00>
<matplotlib.dates.AutoDateFormatter object at 0x000001AC6BF89B20>

这是物理海洋学家 Ken Hughes 在 2021 年发表的博客文章，原文标题为 A better way to code up scientific figures。以 Matplotlib 和 Matlab 为例，强调了模块化思想对于科研作图代码的帮助。我很少看到关于作图方法论的文章，所以翻译出来交流学习。

我画一张出版级别的科研配图一般需要写 100 - 200 行代码，这个长度有点点危险，因为很容易写出能正常运行但又一团糟的东西。如果代码片段都很短还可以从头重写，但如果代码有上千行，提前做好规划会更明智一些。不过在这两种极端情况之间潜藏着另一种吸引人的做法：写出一段当时感觉无比连贯，但以后会让你吃苦头的脚本。

假设你想画一张中等复杂度的图片，类似下面这张：

前言

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 的 pcolor 函数能够绘制由一个个四边形（quadrilateral）单元构成的网格数据的彩色图像，相比绘制等值填色图的 contourf 函数，不会产生过度的平滑效果，能忠实反映像元的数值大小，因而在科学可视化中也很常用。本文并不打算介绍该函数的种种，只想着重讨论网格数据的显示效果、shading 参数发挥的作用，以及 pcolor 和 pcolormesh 这对双胞胎间的差异。本文基于 Matplotlib 3.3.4。

图解网格数据

pcolor 全名 pseudo color，即伪彩色。函数签名为

pcolor([X, Y], C, **kw)

其中 X 和 Y 分别是网格的横纵坐标，C 是网格单元内变量的数值。之所以称之为“伪”，是因为 pcolor 并不像 imshow 那样直接用 RGB(A) 数组表示颜色，而是将 C 的数值归一化之后，在一个颜色查找表中查找对应的颜色，进而用颜色差异表现数值大小（原理详见 Matplotlib 系列：colormap 的设置）。C 数组的形状为 (ny, nx)，X 和 Y 的形状要比 C 大上一圈，即 (ny + 1, nx + 1)，ny 在前表示纵坐标会随数组的行号变动，nx 在后表示横坐标会随数组的列号变动。pcolor 对网格数据的显示效果如下图所示

前言

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