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 包来实现的话，很快就能决定画图思路：

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 坐标系。

简单汇总罗列一下我在网上找到的还不错的 Python 相关资源，包括语言本身以及各种常用库的教程，当然触手可及的官方文档就不收纳了。通通都是免费资源（付费的咱也看不到），分享给有需要的读者。不过互联网资源并非恒久不灭，说不定哪天域名就失效了，或是原作者突然隐藏文章，且看且珍惜吧。

前言

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>

Nicolas Vandeput 发布在 Towards Data Science 上的文章，同时也是其著作《Data Science for Supply Chain Forecasting》中的一章。

为预测任务挑选一个合适的指标并没有想象中那么简单，所以这次我们来研究一下 RMSE、MAE、MAPE 和 Bias 的优缺点。剧透：MAPE 是其中最差的，别用。