简介星载激光雷达CALIPSO-VFM产品数据读取与显示Matlab代码适用于研究地球遥感、气象卫星探测数据的科研人员和学生重点解决垂直特性产品从HDF文件读取到图形化展示的实际问题。压缩包共8个文件由6个Matlab脚本和2个HDF示例数据组成总体仅1.23MB体积轻巧脚本覆盖主程序、图形绘制、动态颜色表创建、类型判断、行块转换等核心环节两个示例数据来自2019年3月不同时刻的标准产品可直接运行查看真实读取与显示结果便于理解数据结构和处理流程。目前已有1501人学习或下载对刚接触星载激光雷达数据处理的人来说是一份容易上手的入门范例。这套代码模块划分清晰读者既能快速复现VFM产品的显示效果也能掌握分块读取、分类着色等关键思路进而扩展针对不同区域或时段的星载激光雷达分析适合在教学、科研或项目预研中使用。轻量体积和模块化设计还可减少调试成本方便在具体课题中直接套用与二次开发。1. 星载激光雷达Calipso-VFM产品最常被问的一句话44字节到底是什么第一次打开 Calipso-VFM 产品的人十有八九会愣一下说好的激光雷达图像读进来不是一张剖面强度图而是一堆 44 字节的掩膜码。这不是数据错了而是 VFMVertical Feature Mask本身就是分类产品——CALIOP 激光雷达把每条视场柱里的目标识别成云、气溶胶、晴空再按类型给出一串编码。做气溶胶和云研究的人用它做类型判据做数据质量控制的人用它剔除无效剖面做算法验证的人拿它当标签真值。这篇文章解决的就是从 HDF4 文件里把这些掩膜读出来、解码、画成能放进报告里的垂直剖面图并把最容易翻车的五个细节一次讲透。2. 先说清楚VFM是什么从L1后向散射到特征掩膜的三个关键概念2.1 VFM在Calipso产品链里的定位不是强度图是分类图Calipso 卫星搭载的 CALIOP 激光雷达向下发射 532nm 和 1064nm 两束激光接收大气后向散射信号。一级产品L1给的是沿轨-高度剖面上的衰减后向散射强度能看到云层和气溶胶层的轮廓但你无法直接从强度值判断“这层是沙尘还是烟”。VFM 属于二级产品它把每一条视场柱、每一个高度 bin 做特征分类输出一个掩膜这个点是不是特征特征类型是云还是气溶胶子类型是什么分类置信度如何。所以读 VFM 时心态要先切换你手上不是一张遥感影像而是一张标注图。每个格点的数值不是物理量而是分类代码。这个区别几乎决定了后面所有处理流程的写法——不能用连续色带渲染不能做平滑插值不能直接把数值当强度用。2.2 文件命名规则与需要读取的HDF4变量VFM 产品以 HDF4 格式分发文件名形如CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf。文件名里值得留意两处一是ValStage1它表示处理阶段不是数据质量等级不能看见 Stage1 就以为是中间产品二是日期时间字段它表示轨段起始的 UTC 时间处理多轨数据时排序别按本地时间排。打开文件后真正要读的变量不多核心是这些变量名作用说明Latitude每个剖面的星下点纬度沿轨方向一维数组Longitude每个剖面的星下点经度与Latitude等长Profile_UTC_Time每个剖面的UTC时间用于横轴和跨轨拼接Lidar_Data_Altitudes高度层网格单位千米从负值到30km附近Feature_Flight_2048_Column_44_Byte特征分类掩膜核心数据44字节逐位编码Feature_Flight_2048_Column_44_Byte这个名字很长拆开看2048 表示沿轨最大剖面数44 表示每个剖面由 44 字节的特征标志组成。有些版本还会出 5km 水平均一的版本变量名里的数字会变成 2060 之类读取逻辑相同。拿到文件后先用下面的脚本把变量名和维度全部打印出来再决定怎么解码。2.3 44字节掩膜的解码思路先拆字节再谈显示44 字节不是 44 个普通整数值而是按位打包的特征标志。每个位或几个连续位对应一个分类结果特征类型、云层类型、气溶胶类型、CAD 分数相关标志、晴空置信度等。直接把这 44 个字节当成 44 个灰度值画图能得到云层轮廓但拿不到类型信息而且不同字节的数值范围差异很大画出来一片噪声。解码的正确思路是先把 HDF4 里的数组按uint8读进来再按位运算提取目标字段。严格来说不同版本V3-01、V3-40、V4-10的位号存在细微调整所以我在实际项目里不会硬编码一套位表而是先把掩膜数组里前几个字节的最小值、最大值打印出来再对照所下载版本的文件说明确认位偏移。这一点后面避坑章节还会展开讲。3. 读取Calipso-VFM一套可直接跑的Python读取流程3.1 安装pyhdf并确认HDF4版本环境准备VFM 是 HDF4 格式这一点是第一个坑。很多新同事拿到.hdf后缀就直接h5py.File(path)立刻报错Unable to open file (file signature not found)。HDF4 和 HDF5 虽然同名文件签名完全不同h5py 只认 HDF5。读取 HDF4 的 Python 库是pyhdf。安装pyhdf在 Linux 下往往要先装 HDF4 运行库。Debian/Ubuntu 系常见做法是先安装系统包再pip install pyhdf。Windows 环境下 pip 直接装 wheel 经常缺依赖我一般建议用 conda 从 conda-forge 安装conda install -c conda-forge pyhdf装完做一次快速验证能成功导入from pyhdf.SD import SD说明环境就绪。这个验证很重要因为pyhdf本身是对底层 HDF4 库的封装库缺失时导入并不会直接报错而是在打开文件时才崩。3.2 最小读取脚本打开文件、列出SD变量、打印shape拿到一个 VFM 文件后我建议第一步永远不要急着读数据而是先列出文件里有哪些科学数据集SD以及每个数据集的维度。这样能避免后面维度顺序搞错导致图横竖颠倒。from pyhdf.SD import SD, SDC path CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf # 以只读方式打开 HDF4 文件 hdf SD(path, SDC.READ) # 列出所有科学数据集及维度信息 for name, info in hdf.datasets().items(): print(name, info) # 关闭文件HDF4 在 pyhdf 里需要显式 end hdf.end()这段代码输出里能看到每个变量的维度元组和数据类型。我见过很多次Feature_Flight_2048_Column_44_Byte的维度在不同版本里发生变化有的版本存成[2048, 545, 44]有的存成[545, 2048, 44]44 这一维是特征标志字节545 是高度层数2048 是沿轨剖面数。如果不先打印维度直接按固定下标取数轻则画图方向错重则整段数据都是错位值。参数说明SD(path, SDC.READ)的SDC.READ是只读模式避免误写原始文件hdf.datasets()返回的是变量名到维度信息字典不是数据本身这一步不会把大数据量装进内存适合快速侦察。3.3 把Feature_Classifier_Flags解码成云/气溶胶/晴空三类掩膜读取核心掩膜并解码我用下面这套流程。代码里先动态把 44 字节那一维挪到最后一维再按位提取特征类型。特征类型的取值按最常见约定0 表示无效或未分类1 表示晴空2 表示云3 表示气溶胶。不同版本可能把 2 和 3 互换但三位分类的思路一致。import numpy as np from pyhdf.SD import SD, SDC path CAL_LID_L2_VFM-ValStage1-V3-40.2018-03-01T00-00-00ZD.hdf hdf SD(path, SDC.READ) # 读取掩膜原始数据 sds hdf.select(Feature_Flight_2048_Column_44_Byte) raw sds.get() # numpy 数组 sds.end() # 统一维序为 (沿轨, 高度, 44字节) arr np.array(raw, dtypenp.uint8) if arr.ndim 3: # 找到长度为44的那个轴并移到最后一维 for i, dim in enumerate(arr.shape): if dim 44: arr np.moveaxis(arr, i, -1) break # 按位提取特征类型用第0字节的最低2位 feature_type arr[..., 0] 0x03 # 构建三分类掩膜0无效, 1晴空, 2云, 3气溶胶 cloud_mask (feature_type 2) aerosol_mask (feature_type 3) clear_mask (feature_type 1) print(原始shape:, raw.shape) print(特征类型取值范围:, np.unique(feature_type)) print(云像元数:, cloud_mask.sum(), 气溶胶像元数:, aerosol_mask.sum()) hdf.end()这里最值得说明的是arr[..., 0] 0x03这一步。arr[..., 0]取出每个格点的第一个字节 0x03是位与操作只保留最低两个二进制位高位的其他特征标志全部置零。为什么用位与而不是直接取整数值因为 44 字节里几乎每个字节都塞了多个标志直接看整数值会被无关高位干扰位与之后才能干净地拿到目标字段。参数说明np.moveaxis(arr, i, -1)的作用是把长度为 44 的轴移动到最后一维这样后面所有下标都用[..., 0]、[..., 1]访问字节位置代码不依赖具体版本里 44 维在哪个位置可移植性好。3.4 数据质量位QA的读取与过滤避免把无效格点画进图里VFM 的 44 字节里除了特征类型还有质量相关标志。我在批量处理时发现如果不过滤质量标志有些夜侧数据或信号衰减严重的剖面会出现一整列“假云”。常见做法是同时读取另一个质量字节把低置信度格点统一改成无效值。质量位的位置同样随版本变化我一般先打印前 8 个字节的统计值观察哪些字节在无效区域恒为 0再确定质量位偏移。下面的代码展示的是典型的过滤流程QA_BYTE和QA_BIT需要按你下载的版本说明调整# 假设质量标志位于第1字节的最高位1表示高质量 QA_BYTE 1 QA_BIT 7 quality_flag (arr[..., QA_BYTE] QA_BIT) 0x01 # 低质量格点与无效格点统一置为0 feature_type[quality_flag 0] 0 feature_type[(arr[..., 0] 0x03) 0] 0 # 重新统计 print(过滤后云像元数:, (feature_type 2).sum())这段代码的核心逻辑是先提质量位再用布尔索引把低质量格点清零。注意不要直接对raw做过滤而是生成新的feature_type数组这样保留原始数据方便回退。参数说明里的两个常量值得单独强调QA_BYTE和QA_BIT在不同 VFM 版本里有差异我习惯把它们放在脚本顶部换版本时只改这两行不碰主体逻辑。4. 显示Calipso-VFM垂直剖面坐标、配色与输出4.1 用matplotlib的pcolormesh画离散色块分类掩膜不能像强度图那样用imshow加连续色带渲染必须用离散颜色映射每个类别固定一个颜色否则云、气溶胶、晴空之间的边界会糊成一片渐变。我用matplotlib的pcolormesh和BoundaryNorm来做因为pcolormesh天然支持非均匀网格而且对离散分类的色块边界控制比imshow干净。import matplotlib.pyplot as plt import matplotlib.colors as mcolors import numpy as np # 沿用上一节的 feature_type, 高度网格和沿轨距离 # heights: 每个高度层中心单位 km # along_dist: 沿轨距离单位 km # 定义类别与颜色0-无效, 1-晴空, 2-云, 3-气溶胶 cmap mcolors.ListedColormap([#f0f0f0, #a6d96a, #1a9641, #d7191c]) # 定义离散边界4个类别需要5个边界值 bounds [-0.5, 0.5, 1.5, 2.5, 3.5] norm mcolors.BoundaryNorm(bounds, cmap.N) fig, ax plt.subplots(figsize(10, 4)) # X轴用沿轨距离Y轴用高度 mesh ax.pcolormesh(along_dist, heights, feature_type.T, cmapcmap, normnorm, shadingnearest) ax.set_xlabel(Along-track distance (km)) ax.set_ylabel(Altitude (km)) ax.set_ylim(-0.5, 20) cbar fig.colorbar(mesh, axax, ticks[0, 1, 2, 3], extendneither) cbar.ax.set_yticklabels([Invalid/Undefined, Clear air, Cloud, Aerosol]) plt.tight_layout() plt.savefig(vfm_profile.png, dpi200)注意feature_type.T这一步feature_type的 shape 是(沿轨数, 高度数)而pcolormesh要求传入的二维数组形状是(高度数, 沿轨数)所以需要转置。这正是前面统一维序的回报——如果开头维度没有理顺这里转置转错一个轴整张图就是镜像的。BoundaryNorm的参数说明bounds列表把数值区间切成 4 段[-0.5, 0.5]对应类型 0[0.5, 1.5]对应类型 1以此类推。这样即使feature_type里有异常数值也不会被归一化拉伸成渐变。4.2 高度轴与沿轨轴的处理让横轴真正代表地面距离VFM 的横轴经常被直接画成剖面序号这在单轨展示时能看但一旦要对比多轨或叠加地理信息序号轴完全不够用。我更习惯把横轴换算成沿轨距离再用副刻度标出经纬度这样读者一眼能看出“这段云层大概在哪个区域”。# 读取经纬度 lat hdf.select(Latitude).get() lon hdf.select(Longitude).get() # 沿轨距离用球面距离累积计算 from numpy import deg2rad, sin, cos, arcsin, sqrt def haversine(lat1, lon1, lat2, lon2): lat1, lon1, lat2, lon2 map(deg2rad, [lat1, lon1, lat2, lon2]) dlat lat2 - lat1 dlon lon2 - lon1 a sin(dlat/2)**2 cos(lat1) * cos(lat2) * sin(dlon/2)**2 return 6371.0 * 2 * arcsin(sqrt(a)) # 计算相邻剖面间的距离并累加 deltas [0.0] for i in range(len(lat)-1): deltas.append(haversine(lat[i], lon[i], lat[i1], lon[i1])) along_dist np.cumsum(deltas)这段代码用 Haversine 公式逐段计算剖面间的地面距离并累计横轴单位是千米。为什么不用np.arange(n_profile)因为卫星在极区附近的星下点轨迹疏密不同剖面序号和实际地面距离不是线性关系尤其在轨道最高纬度的拐弯段序号轴会明显拉伸或压缩距离。用真实距离画出来的云层横向尺度才可信。高度轴的处理相对简单直接用文件里的Lidar_Data_Altitudes就行但要注意这个高度是相对海平面的几何高度单位是千米低层可能从负值开始画图时通常只显示-0.5到20km区间因为 20km 以上只有极其罕见的极地平流层云留太多空白反而削弱图面信息。4.3 配色与图例的工程选择可读性优先不要用彩虹色VFM 显示这块我见过最典型的翻车是直接用viridis或jet渲染掩膜结果云、气溶胶、晴空三类颜色在色带上距离很近打印成黑白稿后完全分不清。离散分类图的核心是类别间的视觉区分度不是美观。我的配色习惯是无效用浅灰晴空用浅绿云用深绿气溶胶用红色。这套方案的逻辑是云和气溶胶是研究者最关注的两类深绿和红在视觉上差异最大晴空作为背景用浅色不会抢前景注意力。如果你需要打印成灰度图建议把云的灰度值压暗、气溶胶提亮确保印刷后仍可分。图例不要用colorbar的默认数值刻度直接改成类别名。上面代码里ticks[0,1,2,3]配set_yticklabels的做法比在图例右上角放一堆数字直观得多也方便直接贴进报告。5. Calipso-VFM读取与显示避坑5条排查记录5.1 现象h5py打开HDF4文件直接报错文件签名找不到原因VFM 是 HDF4 格式h5py 只支持 HDF5两者文件头签名完全不同。这不是文件损坏是库选错了。解决换pyhdf.SD打开。如果项目里已经大量使用 h5py 的读取流程可以把 VFM 先转成 HDF5 或 NetCDF4 再走原流程但我不推荐因为转格式会引入属性丢失和维度重排的额外风险。直接在读取函数里加一个 HDF4 分支更干净。5.2 现象读出来的掩膜数组维度是(545, 2048, 44)按(2048, 545, 44)处理导致剖面方向错乱原因VFM 的维度排列在不同处理版本里不统一代码写死了下标位置换一个文件就出错。解决统一维序后再操作。具体做法是扫描数组各维长度找到等于 44 的那一维挪到末尾再用[..., 0]方式访问字节。这段逻辑我已经放在第 3 章的代码里每次处理新版本文件时先运行维度检查脚本确认后再解码。5.3 现象特征类型提取出来全是0或全是255画出来一片灰原因位偏移用错。VFM 的 44 字节在不同版本里特征类型位所在字节号和位号有调整。V3-01 和 V4-10 的位表就有差异如果直接照抄网上的旧代码很可能取到的是填充字节或质量标志字节。解决先打印掩膜数组每个字节的统计值找到在云区有明显分段特征的字节再用位与提取。我在实际项目里会在脚本顶部维护一个版本到位号的对照字典换数据版本时改一行配置即可。5.4 现象云层剖面画出来是一条紧贴横轴的扁平细线高度结构完全看不出来原因pcolormesh默认的纵横比不是按数据单位等比缩放而显示窗口宽度远大于高度导致垂直方向被严重压缩。解决画出图后先看一眼云层的厚度占比再手动设置ax.set_aspect或用figsize控制画布比例。我一般把figsize设为(12, 4)并将高度轴范围限制在-0.5到20km这样云层的垂直结构能看到又不会因为显示到 30km 高空而浪费大量画布。5.5 现象跨UTC日期的轨段拼接时相邻两轨之间出现一条斜向的“假云墙”原因文件命名里用的是 UTC 日期当地日期跨 UTC 零点时同一天的轨道文件可能分布在两个 UTC 日里按文件名排序会漏掉或错序另外 VFM 每个文件是半个轨道段相邻文件之间本身有重叠或间隔画拼接图时没按时间戳排序就会错位。解决用Profile_UTC_Time作为排序键不依赖文件名日期。拼接前先打印每个文件的起始和结束时间确认顺序连续后再合并拼接图的公共横轴用累计的沿轨距离或直接由时间戳换算。6. 进阶技巧一天多轨段批量出图并输出云/气溶胶占比统计6.1 批处理循环按日期筛选轨道段单轨图画得再好研究里真正需要的是连续一个月甚至一年的统计结果。我处理整日数据时先按轨道段文件名聚合同一 UTC 日期的文件遍历每个文件解码再把同一天的剖面拼到同一个数组里。import glob file_list sorted(glob.glob(CAL_LID_L2_VFM-*.hdf)) daily_masks {} for path in file_list: # 解析UTC日期从文件名第16到25位 date_str path.split(.)[2] # 解码掩膜复用前面封装的函数 feature_type decode_vfm_feature_type(path) if date_str not in daily_masks: daily_masks[date_str] [] daily_masks[date_str].append(feature_type)这段代码的核心是按文件名里的 UTC 日期做分组把一天内的多个轨段收集到列表里。注意这里排序必须用文件名里的 UTC 字段不能用本地时间否则跨 UTC 零点那天会少一段数据。分组之后每个日期对应的列表里每个元素是一轨的二维掩膜下一步就可以逐轨画图或做统计。6.2 统计掩膜占比并输出CSV批处理不能只出图还要把每一轨的云占比、气溶胶占比统计成表格这样后续做纬度带分析时不用重新读数据。import csv with open(vfm_daily_stats.csv, w, newline) as f: writer csv.writer(f) writer.writerow([date, profile_count, cloud_frac, aerosol_frac]) for date_str, masks in daily_masks.items(): # 把同一天所有轨段拼接起来 all_masks np.concatenate(masks, axis0) total all_masks.size cloud_frac (all_masks 2).sum() / total aerosol_frac (all_masks 3).sum() / total writer.writerow([date_str, all_masks.shape[0], round(cloud_frac, 4), round(aerosol_frac, 4)])统计时要注意分母是全部格点数还是有效格点数。我一般用全部格点数同时额外输出一列valid_frac让下游分析能判断统计稳健性——如果某天有效格点占比只有 30%那云占比再高也不具代表性。参数说明np.concatenate(masks, axis0)要求所有轨段的feature_type形状一致包括高度层数。不同版本 VFM 的高度层数可能不同拼之前先统一否则会直接报维度不匹配。6.3 验证方法把统计结果和原始后向散射剖面对照批量出完图后我习惯做一道验证挑一轨数据把 VFM 掩膜叠到同轨的 532nm 衰减后向散射剖面上检查云层和气溶胶层的位置是否与强度高值区吻合。做法是画双面板上面放 L1 强度下面放 VFM 掩膜两图横轴和高度轴完全对齐肉眼核对边界偏移。这一步能快速暴露位号错误、维度错位这类代码问题比直接看统计数字可靠得多。另外一个更廉价的验证是把 VFM 气溶胶出现频率按纬度带统计画成纬度-高度剖面图。真实数据里沙尘层通常出现在北非到中亚的纬度带海洋气溶胶在低纬海洋上空更常见烟羽则常出现在火点纬度附近。如果你的统计图里气溶胶均匀分布在所有纬度大概率是解码位号有问题而不是真实大气如此。这套批处理流程跑过一轮之后我自己的习惯是保留每个版本文件的位号配置、解码脚本和输出图放在同一个目录下换版本时先对照位号再跑批处理避免辛辛苦苦处理完几万条剖面后发现位号用错全部返工。也希望这篇文章能帮你少走这几趟弯路。本文还有配套的精品资源点击获取
