Python 绘制 CALIPSO L2 VFM 产品

灭火器

2023-03-23 更新：

重新简化了代码，详细说了一下原理，详见
https://zhajiman.github.io/post/calipso_vfm/

------------------------------------------------------------
很久以前写过 NCL绘制CALIPSO L2 VFM图像 一贴，现改用 Python 实现，读取并绘制 CALIPSO 反演出的云和气溶胶分类。

数据结构如图所示，代表分类结果的变量 Feature_Classification_Flags 维度为 (3712, 5515)，3712 是沿卫星轨道扫过的点数，5515 是每个点的垂直结构——本来是一个二维的数组，现展平成长为 5515 的一维数组，所以我们在使用时需要恢复它的二维排布。按官网
CALIPSO: Data User's Guide - Data Product Descriptions - Lidar Level 2 5 km Vertical Feature Mask (VFM) Version 3.x Product
的说明，这个二维数组分为垂直三层，每一层的水平分辨率和垂直分辨率都不一样。二维数组沿轨道方向长 5 km，高 30 km，文件中自带的经纬度和时间对应 5 km 中点处的值。实际使用时为了方便起见，取每一层的第一条廓线来代表此层的结果，也就是图中标红的那三条。于是每个点的二维数组简化为垂直分辨率不一的一条廓线。

廓线还不能直接画，每个数组元素的数据类型是 ushort，在 numpy 里对应于 np.uint16，也就是一个 16 位的无符号整数。以二进制形式显示一个元素的值，如图所示从最右往左分为 7 组，CALIPSO 将分类结果编码到了这 7 组比特中。例如最右边第一组代表分类的主要结果，如果其值为 011，也就是十进制的 3，那么就代表气溶胶。Python 中可以通过位运算提取这些信息，具体解读详见官网。
这里使用 pyhdf 包读取 CALIPSO 的 HDF4 文件，代码中用到的 map_funcs 请见 Cartopy 系列：从入门到放弃 一文最后一节。

使用 MATLAB 或 IDL 的同学可以直接调用官网提供的代码。
参考链接：
https://hdfeos.org/zoo/index_openLaRC_Examples.php
http://www.meteothink.org/examples/meteoinfolab/satellite/calipso.html
下面以 2021-03-15 中国北方发生的强沙尘暴为例，VFM 垂直剖面图中可见大量气溶胶分布。
为了避免用 add_axes 方法调来调去，这里用 Axes 的 get_position 方法得到地图的实际大小后，定量控制各个 Axes 间的位置（代码里相关的部分还可以参考 添加子图的便利函数 一帖进行简化，这里就懒得改了）。
利用线性插值的方法得到剖面图中各个刻度的经纬度及其标签。

mico2008

{:5_213:}{:5_213:}{:5_213:}{:5_213:}

独孤酒见

楼主太牛了。
膜拜一下。
楼主再辛苦完善下把x轴的位置标注为下面的这种形式呗
lat     lat
lon    lon

灭火器

独孤酒见 发表于 2021-12-18 09:53
楼主太牛了。
膜拜一下。
楼主再辛苦完善下把x轴的位置标注为下面的这种形式呗

已编辑（      ）

931404656

发错了不能撤销

独孤酒见

灭火器 发表于 2021-12-18 10:10
已编辑（      ）

牛人。
太感谢了。
辛苦了。
周末快乐。

uranus

太牛了，这个数据存的也是奇怪，5515的划分真是清奇

Liujiaxin

学习到了，多谢

white_lcz

谢谢楼主！学习很多，获益匪浅！
请问楼主有没有想过使用其中的quality flag...来筛选出高质量的数据？
感谢，顺祝新年快乐！

灭火器

white_lcz 发表于 2022-2-7 11:56
谢谢楼主！学习很多，获益匪浅！
请问楼主有没有想过使用其中的quality flag...来筛选出高质量的数据？
...

质量控制信息很容易提取

1. f = ReaderCAL(filepath)
   
2. ftype = f.read_ftype()
   
3. ft = ftype[:, :, 0]
   
4. qa = ftype[:, :, 1]

复制代码

其中 qa 即 ft 的质量控制标记。不过我还没有在研究中用到，所以对质量相关的信息没什么了解。