R语言处理CALIPSO气溶胶廓线：HDF4读取与剖面图绘制指南

大概两年前我第一次拿到CALIPSO L2 5km气溶胶廓线产品想用R把532nm消光系数画成沿轨剖面图。当时的第一反应是这不是跟处理普通遥感栅格差不多嘛读进来画个热力图半小时搞定。结果光是让R正确读出一个HDF4文件就折腾了大半天后面还有维度方向、质量控制、云层污染这一串坑等着。这篇文章就把我后来沉淀下来的完整处理流程写出来包括读取方式选型、字段地图、质量控制逻辑、剖面图设计以及几个最容易让人怀疑人生的bug排查思路。适合正在用R处理CALIPSO、CALIOP或者其他激光雷达廓线产品的朋友尤其是第一次接触HDF4格式、还没被消光系数折磨过的人。1. 先搞清楚你手里到底是什么L2产品与532nm消光系数1.1 消光系数不是直接测出来的是反演出来的532nm消光系数是CALIPSO二级产品里最常被使用的核心变量之一但很多初学者把它当成激光雷达直接测量的量这是个容易造成误解的地方。激光雷达直接测量到的是衰减后向散射信号也就是L1产品里的Attenuated Backscatter它是光和大气相互作用之后的综合结果。而L2产品里的消光系数是在假设激光雷达比Lidar Ratio的前提下通过约束反演把后向散射信号转换成消光廓线。消光系数单位是km⁻¹物理意义是光在单位距离上被衰减的比例数值越大表示该高度层的气溶胶或者云越浓。532nm这个波段在CALIPSO里主要用来反演气溶胶光学特性因为短波对细粒子更敏感。也正因如此这个变量做出来的图纵向看能反映气溶胶的垂直分布结构横向看能反映沿卫星轨道的空间变化是一张信息量很大的图。理解这一点很重要因为后面做图时你会遇到消光系数出现负值某个廓线完全不可靠等情况。如果只是把数据读出来填进矩阵再画图很容易被反演噪声带偏。数据文档里通常会给有效取值范围比如0附近的小正值才是可信的明显偏离物理量级的负值或者超大正值都要当异常处理。1.2 文件是HDF4容器R默认读不了CALIPSO L2文件是老牌的HDF4格式。这里说的HDF4和R语言里常见的HDF5并不是一回事很多人在这一步就卡住了。ncdf4处理的是NetCDFraster虽然支持部分HDF但主要是GDAL驱动下的HDF4/5子数据集面对CALIPSO这种非地理网格化的轨道廓线数据经常表现得非常别扭。HDF4这种格式本质上是一个容器里面可以塞多个数据集SDSScientific Data Set每个数据集都有自己的名字、维度、属性。CALIPSO L2文件里消光系数、质量控制标记、经纬度、时间、特征类型等变量都是独立的SDS没有统一的地理参考坐标系因为卫星沿轨扫描并不是一个规则网格。所以处理它的关键思路是不要试图把整个文件当作一张遥感影像来读而是要按变量取数。1.3 处理前先画一张字段地图我强烈不建议拿到文件就直接写R代码。先花十分钟把文件结构摸清楚后面能省几个小时。用工具查看文件内容之后你需要确认下面这些字段是否存在、命名是什么、维度是什么字段类型常见命名作用消光系数Extinction_Profile_532核心数据km⁻¹质量控制Extinction_QC_Flag判断每个点是否可信特征类型Feature_Type区分云、气溶胶、无效层云气溶胶判别CAD_Score云/气溶胶判别得分高度Height / Initial_Height垂直层坐标时间Profile_UTC_Time每个廓线的时间经纬度Latitude / Longitude每个廓线的水平位置不同版本号的文件字段命名可能有细微差别比如有的版本里会有Total_Number_Height_Levels定义层数有的版本把高度信息放在数据集的属性里而不是单独的数据集里。所以无论你看谁的博客、用谁的代码第一件事永远是用工具列出实际字段而不是把网上抄来的路径硬套进去。这一步我称之为字段地图有了它后面的读取代码就是填空而已。2. 三种读取路径对比别一上来就死磕HDF42.1 首选路径rhdf5直接读取多数情况能开我要先说一个反直觉的事rhdf5这个包虽然名字里写的是HDF5但在很多预编译环境里它依赖的底层库带了HDF4兼容支持。也就是说你完全可以先用rhdf5试试直接打开CALIPSO的HDF4文件。我试过的多台机器上大部分L2文件都能被它正常读取只有少数版本或损坏文件会报错。library(rhdf5) f - CAL_LID_L2_05kmAPro_V4-20_2021-001T00-00-00ZD.hdf # 查看文件内部结构 h5ls(f) # 读取核心数据 ext - h5read(f, /Extinction_Profile_532) qc - h5read(f, /Extinction_QC_Flag) lat - h5read(f, /Latitude) lon - h5read(f, /Longitude) time - h5read(f, /Profile_UTC_Time) # 查看维度 dim(ext) length(lat)文件结构里通常会看到很多/开头的路径这就是SDS在HDF5视图下的位置。读取之后发现维度和你预期不符先别慌大多数情况下是行和列的顺序问题后文有处理办法。rhdf5没有正式公开支持HDF4所以如果在你的环境里打开失败、报类似unable to open file的错误不要死磕跳到第二条路径——用GDAL做转译。这并不是你的代码有错而是你的底层HDF5库在编译时没有启用HDF4支持。2.2 稳妥后备GDAL命令把子数据集拆成R能吃的格式GDAL工具链处理HDF4的能力远比R自带生态成熟。操作分两步先在命令行里用gdalinfo查看文件内部有哪些子数据集gdalinfo CAL_LID_L2_05kmAPro_V4-20_2021-001T00-00-00ZD.hdf输出里会列出类似HDF4:filename.hdf:/Extinction_Profile_532这样的子数据集路径。确认路径名之后用gdal_translate把关心的数据集单独抽出来转成NetCDF或GeoTIFFgdal_translate -of netCDF \ HDF4:CAL_LID_L2_05kmAPro_V4-20_2021-001T00-00-00ZD.hdf:/Extinction_Profile_532 \ extinction_532.nc转出来的NetCDF可以直接用ncdf4读取维度关系一目了然。如果连gdal_translate都不想用还可以用Panoply这类图形工具手动导出CSV但只适合单条轨道、一次性出图批量处理不推荐因为效率太低。2.3 为什么我基本不用raster直接读有人可能会问raster包不是也能读HDF吗是的raster能读但它隐含着一条规则把子数据集当作风靡一时地理栅格行和列会被理解为经纬度网格。CALIPSO的消光廓线矩阵本质上是高度层 x 沿轨廓线的关系不是经纬度规则网格。一旦让raster介入它可能强行赋予地理参考导致后面坐标全乱。我用过一次之后就再也不这么干了。归根结底读取方式的选型逻辑很简单优先用能保持数据原始组织方式的接口rhdf5或ncdf4尽量避免让栅格类专业工具二次解释数据。3. 数据清洗与坐标重建出图前最脏最累的环节3.1 维度顺序一定要用它自己的定义来核对我先说一个出现频率极高的问题维度反了。HDF文件里的数据矩阵不同变量组织方式不一样不能想当然地认为dim(ext)的第一个维度就是高度。我的习惯是读入之后立刻打印维度并和文件里记录的维度名称做比对。h5ls(f)输出里通常会写明每个数据集的维度名如果消光系数字段标识为[height, profile]那就说明第一维是高度、第二维是廓线。如果标识反了后面所有处理逻辑都要跟着翻转。这里给一个典型的重排代码把高度 x 廓线变成R里处理起来更顺手的形式nlev - dim(ext)[1] nprof - dim(ext)[2] # 假设第一维是高度第二维是廓线 ext_by_prof - t(ext) qc_by_prof - t(qc) # 转成 tidy 格式 df - data.frame( profile rep(1:nprof, each nlev), height rep(1:nlev, times nprof), ext as.vector(ext_by_prof), qc as.vector(qc_by_prof) )如果你的数据维度说明是反过来的把t()去掉即可。关键规律是先确认再重排不要盲目相信网上的任何代码。3.2 无效值和质量控制分三步过滤直接读出来的矩阵里会混着各种无效点、可疑反演值。我习惯按下面三步依次过滤。第一步处理填充值。HDF里的无效点通常用一个很大的负值表示比如-9999。先把它替换成NAdf$ext[df$ext 0] - NA第二步应用质量控制标记。Extinction_QC_Flag是一个整数标记不同数值段代表不同质量等级。文档里通常会写清楚哪几个值是完全有效“部分可用”哪几个是不可靠。我常用的做法是只保留质量标记为完全有效的点# 以常见约定为例0 通常表示高质量 df$ext[df$qc ! 0] - NA第三步也是很多人容易忽略的一步用Feature_Type做语义过滤。消光廓线里不是所有层的值都代表气溶胶云层同样会产生强大的消光信号。如果你只关心气溶胶那就要把云层剔除。Feature_Type里常见编码是1代表云、2代表气溶胶、3代表平流层气溶胶但这些编码在不同产品版本里可能不同。用文档核对之后构建掩膜feature_type - h5read(f, /Feature_Type) # 保留气溶胶相关层 keep_mask - feature_type %in% c(2, 3) df$ext[!keep_mask] - NA这三步做完数据才算干净后面画出来的图才有解释价值。3.3 高度不是序列号是用Initial_Height换算出来的很多第一次处理CALIPSO的人会把矩阵行号当成高度层序号来画图画出来的纵轴是整数层号别人根本没法读。正确做法是拿到每个高度层的实际物理高度。CALIPSO L2产品的高度信息有两种存在形式一种是单独的Height数据集另一种是放在Extinction_Profile_532属性里的Initial_Height。如果是后者每个廓线的起始高度可能相同也可能随轨道渐变比如地形起伏导致信号起点变化需要用h5readAttributes把属性取出来attrs - h5readAttributes(f, /Extinction_Profile_532) initial_height - attrs$Initial_Height # 单位通常是km layer_width - attrs$Width_of_Layer # 每层厚度通常0.03km或0.06km height_vec - initial_height (0:(nlev - 1)) * layer_width然后把height_vec按廓线重复拼进长表数据框。这里特别注意一定要确认高度单位是公里还是米文档或属性里有写。我见过有人把单位弄错画出来的剖面图高度轴从0到数百公里室内汇报直接被问住。时间坐标同理。Profile_UTC_Time字段通常是一个相对某个参考时间点的时间偏移量可能是秒也可能是毫秒。读入后用起始时间加偏移量转成标准时间格式# 把整型偏移按文档说明转成秒 start_time - as.POSIXct(1993-01-01 00:00:00, tz UTC) df$datetime - start_time as.numeric(time)[df$profile]这里的时间偏移基准要看文档说明不同产品可能不同不要照抄任何博客的固定数值。这个坑我踩过之后总结了一句话数据里凡是带单位或者基准说明的东西都要当成假设去验证而不是当成事实去使用。4. 从矩阵到剖面图可视化方案的取舍4.1 剖面图的设计逻辑对数色标与高度裁剪数据清洗完之后绘图本身并不复杂但要注意几个视觉设计上的关键点。第一一定要用对数色标。消光系数的数量级跨越非常大洁净大气层可能只有0.0001 km⁻¹而云层或者强沙尘层可以到0.1甚至1以上。线性色标下大气边界层的气溶胶会被云层或者个别强信号直接压成一片蓝色什么都看不出来。取对数之后低值区的丰富结构才能显示出来。第二要对色标范围做截断。log10之后范围通常集中在-4到-1之间个别强信号或者未过滤掉的噪声可能到0以上。如果不截断色标两端会被极端值占据大部分数据只能挤在中间很小的色带里。我会用limits加oob参数配合让超出范围的值用边缘色表示。第三高度范围要裁剪。5km分辨率消光廓线的垂直范围能到30km左右但平流层中高层的信噪比很差画出来通常是大片噪声。我一般只保留0到15km这个范围内气溶胶和云的垂直分布信息已经比较完整图面也更聚焦。4.2 ggplot2实现一张可发布的剖面图下面这版代码是我目前用得最顺的模板包含对数色标、高度裁剪、云层叠加library(ggplot2) library(dplyr) plot_ext_profile - function(df, cloud_df NULL) { p - df %% filter(height 15, ext 0, is.finite(log10(ext))) %% ggplot(aes(x datetime, y height, fill log10(ext))) geom_tile(width 2) scale_fill_gradientn( colours c(#3B4CC0, #2E8CBE, #68B6A8, #E0DF9A, #F0B32A, #E64B35), limits c(-4, -0.5), oob scales::squish, name log10(ext/km-1) ) scale_y_continuous(limits c(0, 15), expand c(0, 0)) labs(x UTC time, y Altitude (km)) theme_minimal(base_size 14) if (!is.null(cloud_df)) { p - p geom_point( data cloud_df, aes(x datetime, y height), colour white, size 0.2, shape 4 ) } p }geom_tile的width参数我设成了2秒或对应的轨道距离避免瓦片之间出现细缝。云层叠加用的是Feature_Type 1的层用白色叉号标出来这样读者能看到气溶胶信号周围云的位置解释数据时更有底气。4.3 批量处理多轨数据函数化是关键做区域分析的时候很少只处理一条轨道有时候一个月的数据就要几十个文件。我强烈建议把读取、清洗、绘图封装成函数再用lapply或purrr循环处理。文件名里通常嵌着日期和时间信息写个小解析函数提出来library(purrr) library(stringr) process_one_file - function(f) { ext - h5read(f, /Extinction_Profile_532) # ... 中间清洗步骤与上文学到的坐标重建 ... df } file_list - list.files(data/calipso/, pattern *.hdf$, full.names TRUE) all_df - map_dfr(file_list, process_one_file, .id track_id)合并后的数据框加上track_id列做分面图ggplot(all_df, aes(x datetime, y height, fill log10(ext))) geom_tile() facet_wrap(~track_id, ncol 1, scales free_x) scale_fill_gradientn( colours c(#3B4CC0, #2E8CBE, #68B6A8, #E0DF9A, #F0B32A, #E64B35), limits c(-4, -0.5), oob scales::squish )批量出图的成果放到汇报里就是很直观的区域气溶胶垂直结构时间序列图。4.4 出图前检查清单我给自己定了一个出图前检查习惯每张图发布前都过一遍纵轴单位是km且高度范围合理横轴是真实时间而不是廓线序号log10之前所有非正值都被过滤或赋值为NA色标范围不是自动默认的而是截断过的云层和气溶胶至少看一眼特征类型分布确认语义正确这套清单看起来琐碎但真的能拦住很多数据看起来没问题但图一发表就被审稿人问住的情况。5. 高频踩坑的完整排查链路从全灰图到锯齿剖面5.1 读取出来全是灰色空图先查子数据集路径再查维度现象费劲读出来的矩阵画图后全图一个颜色或者95%以上是空白。排查链路是这样走的第一步查看sum(is.na(df$ext))占总数比例。如果占比极高说明大多数值本来就被判为无效问题大概率出在质量控制过滤太激进或者你读的根本不是消光系数数据本身。第二步回到h5ls(f)输出确认读取路径是对应Extinction_Profile_532而不是别的变量。这种情况我经历过一次读取时路径少了一层结果把某个属性数组当成主数据读了出来矩阵维度和坐标对不上画出来自然全是乱的。第三步检查QC过滤是否合理。如果你用的质量标记约定和当前产品版本不一致可能把有效数据全过滤掉了。这时候干脆先不完全过滤只处理填充值画一张裸图看看数据本身是否有信号再逐步加回过滤条件。5.2 剖面图高度轴出现倒置或锯齿现象剖面图显示出的气溶胶层是上下颠倒的或者某些高度层出现横向锯齿状条纹。先说倒置。CALIPSO的廓线层有的是从低海拔到高海拔排列有的在高海拔产品里是从高到低排列。如果你发现图里信号最强的层出现在高空而不是边界层极有可能就是高度向量方向反了。处理方法很简单if (height_vec[1] height_vec[2]) { height_vec - rev(height_vec) # 矩阵每一列也要反转 ext_by_prof - apply(ext_by_prof, 2, rev) }锯齿条纹更多见的原因是高度坐标重建时没有处理好不同廓线的起始高度。如果每个廓线的Initial_Height在变化而你硬用了同一个高度向量去映射所有廓线那么同一气溶胶层在不同廓线里会被画到不同高度视觉上就出现锯齿状。解决办法是构造完整的高度矩阵而不是单一向量height_mat - matrix( rep(initial_height_vec, each nlev) rep((0:(nlev - 1)) * layer_width, times nprof), nrow nlev, ncol nprof )这部分数据量会变大但这是正确表达轨道数据结构的方式。5.3 云层信号压制气溶胶色标现象整张图色标范围完全被几个红色亮点控制气溶胶信号在上百公里宽的图面上只占很小一片蓝绿色区。这就是典型的云层污染。云层的消光系数往往比气溶胶高一两个数量级如果不过滤掉对数色标的动态范围也会被拉扯。解决思路不是调色标而是回到数据语义层只保留气溶胶层或者把云层单独作为叠加标记显示。我用的是后者因为做大气研究时云层位置本身也是有价值的信息。如果必须把云和气溶胶画在一起那就单独用另一个色标编码云层或者把云层画成透明度较高的白色覆盖层而不是让它们在一个连续色标里互相挤压。5.4 白天轨道噪声大什么时候应该直接放弃现象白天轨道的剖面图里整条图都布满细碎的高值亮点尤其是低层看起来像撒了一把胡椒面。这是激光雷达白天信噪比差导致的太阳背景光会淹没微弱的后向散射信号。CALIPSO白天轨道的数据质量在低层尤其糟糕。我见过不少人花大量时间调参数试图把白天轨道的图修好效果都很有限。更务实的做法是做气候或月平均分析时尽量只用夜间轨道或者在合图中对白天轨道单独标注。有一回我做某区域三月份的气溶胶垂直结构图混用了白天和夜间轨道图面看起来信息量很大但细看之后发现很多结构其实是噪声差点把结论带偏。后来我做了个简单实验同一天白天轨道和夜间轨道的消光剖面放在一起对比信噪比差距一目了然从此之后我处理L2消光系数时默认夜间轨道优先。最后再分享一个我个人的习惯处理CALIPSO这类格式老、结构怪的数据我最怕的不是代码写不出来而是代码跑通了但结果其实是错的。所以我每写一步处理逻辑都会顺手用最基础的方式验证一下——比如读出一个廓线的消光值和Panoply里同一位置肉眼对比一下数量级比如画图之前先打印高度向量前三个值确认是0.0、0.03、0.06而不是倒序。听起来很笨但正是这些笨功夫让我省下了无数在汇报桌上被质疑的时间。R处理消光系数的路子不算多但每一步都走扎实的话画出来的图足以支持你从现象观察到结论输出。