以前处理功率时间序列时经常遇到一大段时间里功率值虽然没有缺失，但全是零的异常情况，为了找出这些连续为零的时段，当时设计了一个统计序列里零值连续出现次数的函数，效果如下图所示：
输入序列是
series = np.array([0, 0, 1, 2, 1, 0, 0, 0, 0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 0, 0, 3, 4, 3, 0])
其中有四段零值，长度依次为 2、4、6、1。输出序列与输入序列等长，输入序列中非零位置的数值为零，零值位置数值为零值连续出现的次数。
最早笔者用 Python 画中国地图时，会准备 bou2_4p.shp 文件，然后封装一个读取 shapefile 并添加到 GeoAxes 上的函数，别的项目要用时就把数据和函数复制粘贴过去。Cartopy 系列：从入门到放弃 里就是这么做的。
后来工作中用到了 Clarmy 开发的 cnmaps 包，只用两行就能快速绘制地图，非常方便。同时萌生了自己实现一个功能类似的包的想法，遂开发出了 frykit。
ERA5 的 NetCDF 文件或卫星的 HDF 文件为了压缩文件体积会用 16 位整数存储变量，读取时跟属性里的
add_offset和scale_factor做运算恢复成 64 位浮点数。如果你是用 Python 的 NetCDF4 或 xarray 包处理 NetCDF 文件，甚至都不用关心这些细节，它们默认会帮你解包成浮点数。问题是，如果自己也想用这种方法压缩数据，那么add_offset和scale_factor该如何设置，压缩率能有多高，又会损失多少精度呢？一番搜索后发现 Unidata Developer’s Blog 上的博文 Compression by Scaling and Offfset（原文标题确实把 offset 拼错了）清晰地介绍了压缩的原理和参数选择，现翻译前半部分，后半部分关于 GRIB 压缩的看不懂感觉也用不上，偷懒不翻了。
今天来深入了解一下存储浮点数据时如何指定所需的精度，抛弃那些对于精度来说多余的比特。这些多余的比特往往很随机所以不可压缩，导致标准压缩算法的效果有限。需要注意这种操作是一种有损压缩。
人眼可见色域在色度图中表现为彩色的马蹄形，单色光（monochromatic light）的颜色对应于马蹄的弧形边界。本文想将单色光的颜色按波长线性增大的顺序一字排开，用类似彩虹渐变图的形式展示单色光光谱。用 Python 的 Matplotlib 包来实现的话，很快就能决定画图思路：
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 坐标系。
简单汇总罗列一下我在网上找到的还不错的 Python 相关资源，包括语言本身以及各种常用库的教程，当然触手可及的官方文档就不收纳了。通通都是免费资源（付费的咱也看不到），分享给有需要的读者。不过互联网资源并非恒久不灭，说不定哪天域名就失效了，或是原作者突然隐藏文章，且看且珍惜吧。
CALIPSO 卫星的 L2 VFM（Vertical Feature Mask）产品根据激光的后向散射和消光信息，将激光通过的各高度层分类为云或气溶胶。该产品在现实中的表现如下图所示：卫星一边在轨道上移动一边向地面发射激光脉冲，相当于在地面上缓缓拉开一幅“画卷”，VFM 描述了“画卷”上云和气溶胶的分布和分类情况。
处理 VFM 产品的难点在于：
(N, 5515) 的形状，N 表示卫星移动时产生了 N 次观测，但 5515 并非表示有 5515 层高度，而是三种水平和垂直分辨率都不同的数据摊平成了长 5515 的数组。因此处理数据时需要参照文档的说明对 5515 进行变形。网上能找到的代码有：
笔者也曾写过两次教程：
本文是对旧教程的翻新，会对 VFM 数据的结构进行更多解释，对代码也进行了更新。本文使用 pyhdf 读取 HDF4 文件，用 Matplotlib 3.6.2 画图。为了方便画图，用了一些自制的函数（frykit）。虽然基于 Python，但希望能给使用其它语言的读者提供一点思路。
完整代码已放入仓库 calipso-vfm-visualization。