Matlab批量读取GOCI卫星.he5文件并提取TSS浓度数据的方法

如果你手头刚好有一批GOCI卫星的水色产品文件后缀是.he5需要一次性读进来、提取里面的总悬浮物质浓度TSS数据再按时间顺序画成变化过程图那这篇内容应该可以帮你省不少摸索时间。我用Matlab处理GOCI数据有几年了从最初单文件手动调到后来批量跑几十个时相中间踩过的坑基本集中在文件结构不熟、批量组织混乱、比例因子忘乘这三大类。这篇文章把整套流程拆开讲清楚包括如何用h5info定位TSS数据、如何写批量循环、如何做掩膜和配色最后也会给一份可以直接改路径拿去跑的Matlab代码。适合刚接触水色遥感的同学也适合已经被.he5文件结构折磨了一阵子的从业者。1. 动手之前先搞清楚GOCI和TSS是什么1.1 GOCI数据到底特别在哪GOCI是搭载在通信、海洋、气象卫星上的静止轨道海洋水色成像仪它跟极轨卫星最大的区别是能“盯着一个地方看”。极轨卫星一天过境一两次而GOCI可以在白天对同一海域进行多次观测时间分辨率能到小时级。这意味着什么呢如果你做的是潮周期内的悬浮泥沙变化、突发污染扩散、或者近岸水体短期演变分析GOCI的时间密度是很大的优势。TSSTotal Suspended Solids即总悬浮物质浓度反映的是水体内悬浮颗粒物的含量单位通常是g/m³或mg/L。近岸水质评价、河口泥沙输运、赤潮跟踪TSS都是一个硬指标。我所接触的GOCI L2产品里TSS并不是原始遥感反射比而是已经经过大气校正和反演算法处理的二级产品。你拿到文件之后不需要自己从头算只需要把产品变量读取出来处理好比例因子和无效值然后做空间展示或时间序列分析就好。当然前提是你得知道怎么把这个.he5文件读进来这就是接下来的重头戏。1.2 TSS产品的数据链路GOCI数据处理系统GDPS是官方提供的数据处理平台用户从数据共享平台拿到的一般是L1原始数据L1经过GDPS处理之后输出L2产品其中就包含TSS、叶绿素a浓度、有色溶解有机物、总吸收系数等若干波段。输出文件通常以HDF5格式存储而GDPS有个习惯文件后缀用的是.he5而不是.h5。这个细节很关键因为很多人在文件浏览器里看到.he5第一反应是“这什么格式”其实它就是HDF5Matlab的h5系列函数完全可以直接读。以我实际接触过的产品为例文件内部结构通常走HDF-EOS规范。所谓HDF-EOS可以理解为NASA对HDF5做了一套面向地球观测的结构约定把数据组织成“Swath”或“Grid”Swath内部再分成数据场、地理定位场、属性等节点。打开之后你经常能看到类似这样的路径/HDFEOS/SWATH/GOCI/Data Fields/TSS/HDFEOS/SWATH/GOCI/Data Fields/CHL/HDFEOS/SWATH/GOCI/Geolocation Fields/Longitude/HDFEOS/SWATH/GOCI/Geolocation Fields/Latitude并不是所有产品都完全一致不同处理版本内部路径可能略有差别所以最稳妥的做法不是凭记忆硬编码路径而是先用h5info看一遍文件结构。我在后面会专门演示这一步。了解清楚这个结构后续批量处理才不心虚。2. 环境准备与批量读取方案2.1 Matlab读HDF5的工具与准备Matlab对HDF5的支持相当完善不需要额外装气象或海洋工具包只需要基础版就能调用h5info、h5read、h5disp这几个核心函数。h5info用来查看文件结构返回一个结构体里面包含Groups、Datasets、Attributes等字段h5read用来按路径读取数据h5disp则是把文件结构直接打印到命令行适合快速扫一眼。有一点需要提醒旧版本Matlab比如2015年前后对HDF5的某些压缩算法支持不是很好如果读取写到一半报错提示“filter”相关优先考虑升级Matlab版本或改用其他读取工具。另外文件名路径里不要带中文和特殊符号个别版本在Windows下对中文路径敏感读不出来的时候先把路径改成纯英文试一次。除了基础环境我建议在脚本开头固定设置这样一段clear; clc; close all; warning(off, all); % 关闭非关键警告但不建议全程关闭批处理时warning太多会严重拖慢速度我一般在调试阶段不开这条跑批量的时候才关掉。2.2 批量读取.he5文件的两个核心套路批量读取的核心就是把目录下所有.he5文件找出来再逐个循环。最常用的方法有两种第一种是用dir加通配符第二种是用文件夹筛选函数。实测算起来dir配合fullfile最稳跨电脑也不会因为斜杠方向出问题。下面这段代码就是批量读取的骨架你可以直接当成模板fileDir D:/GOCI/TSS/20230731; fileList dir(fullfile(fileDir, *.he5)); fileList fileList(~[fileList.isdir]); % 排除文件夹 fprintf(共找到 %d 个文件\n, numel(fileList)); for k 1:numel(fileList) fileName fullfile(fileDir, fileList(k).name); fprintf([%02d] %s\n, k, fileList(k).name); % 这里写你的读取与处理逻辑后面会补全 % tss h5read(fileName, tssPath); end这里有个容易被忽略的点dir返回的结构体里第一个元素有时候是当前目录虚拟项isdir为true如果你用通配符’*.he5’这种情况通常不会出现但加上排除文件夹那行可以防御性避免出错。2.3 文件名里的时间信息怎么利用GOCI的文件名一般长这样例如“COSI_GOCI_L2_20160731024700.he5”中间那串数字就是观测的UTC时间年、月、日、时、分、秒基本都齐了。画TSS变化过程最关键的一步就是把这个时间解析出来并用它来给数据排序和作为横轴。我自己常用的解析方法是正则加datestrtokens regexp(fileList(k).name, \d{14}, match); if ~isempty(tokens) timeStr tokens{1}; timeNum datenum(timeStr, yyyymmddHHMMSS); timeVec(k) timeNum; end注意GOCI文件里记录的时间是UTC如果海域位于东八区实际当地时间需要加8小时。要不要当地时取决于你的应用场景但如果做潮汐相关的分析时间基准必须统一否则相位会差八小时结论可能完全相反。还有个小技巧预处理阶段就把所有文件路径和时间戳写进一个表格结构体后面读数据时按时间排序。这个习惯能让你在处理几十上百个文件时思路清晰meta struct(name, {}, path, {}, timeNum, {}); for k 1:numel(fileList) ... meta(k).name fileList(k).name; meta(k).path fullfile(fileDir, fileList(k).name); meta(k).timeNum datenum(tokens{1}, yyyymmddHHMMSS); end [~, idx] sort([meta.timeNum]); meta meta(idx);3. 提取TSS数据并做预处理3.1 找到TSS数据的真实路径这一步是整个流程里最容易出岔子的地方。一个优质的TSS产品文件里除了TSS之外通常还有叶绿素、CDOM、云掩膜、质量标记等一大堆字段路径结构也分好几层。你要是靠回忆去找往往看了半天也没找到。我的习惯是拿到一个陌生文件后第一件事就是跑h5disp把结构打印出来然后从输出里直接搜TSS关键字。比如这样info h5info(D:/GOCI/TSS/20230731/COSI_GOCI_L2_20230731224700.he5); disp(info);然后命令行里会显示出整个树形结构。找到类似‘/HDFEOS/SWATH/GOCI/Data Fields/TSS’这样的路径后我建议把路径写成一个变量放在脚本顶部不要散落在各处硬编码tssPath /HDFEOS/SWATH/GOCI/Data Fields/TSS; lonPath /HDFEOS/SWATH/GOCI/Geolocation Fields/Longitude; latPath /HDFEOS/SWATH/GOCI/Geolocation Fields/Latitude;如果你用我上面提到的方法查看多个文件发现不同产品的路径有差异最快的办法是先写一个通用函数自动扫描文件内所有含TSS的数据集路径。虽然这次你可能不需要但等到你换一个版本的数据源时能省去重复人工排查的功夫。3.2 比例因子、有效范围和FillValue处理HDF5产品里的数据为了压缩体积存储时通常不是浮点数而是int16或uint16的整型。这意味着直接读出来的数值并不是物理量本身需要根据文件属性里的ScaleFactor和Offset做转换。GDPS产品常见的比例因子可能是0.1或0.01具体值取决于变量和版本千万别想当然。我调试时遇到过一个很典型的情况读出来的TSS矩阵里数值从0到60000多随手一画图就黑漆漆一片后来仔细翻了属性才看到单位是“0.1 g/m³”也就是说实际浓度是存储值乘以0.1。这一步没做对后续所有分析全部白搭。推荐的转换代码如下% 读取原始整型数据 tss_raw h5read(fileName, tssPath); % 读取属性里的比例因子和偏移量 scale h5readatt(fileName, tssPath, scale_factor); offset h5readatt(fileName, tssPath, offset); fillValue h5readatt(fileName, tssPath, _FillValue); % 转成浮点并应用比例因子与偏移量 tss_double double(tss_raw) * scale offset; % 把填充值和负值统一替换为NaN tss_double(tss_raw fillValue | tss_double 0) NaN;属性名不一定百分百固定有的文件里写的是ScaleFactor有的写的是Slope有的干脆没有offset。稳妥做法是用try-catch读属性读不到就默认scale为1、offset为0。这个处理逻辑在批量循环里特别重要因为个别文件可能属性缺失不做保护整个循环就会中断。3.3 陆地、云与无效像元的掩膜提取出有效浓度后还需要把不想显示的像素盖掉。第一类是陆地像元水色产品的范围往往已经裁到海区但近岸处仍可能包含陆地第二类是云像元GOCI产品里通常有质量标记或云掩膜变量第三类就是前面说的填充值。如果不做掩膜图上一片陆地反射值会被当成高浓度水体配色比例尺会被拉得很离谱。如果你能找到云掩膜变量读取后把非水体像元设为NaNmask h5read(fileName, maskPath); validMask (mask 1); % 具体有效编码需看属性 tss_double(~validMask) NaN;如果产品里没有云掩膜也可以用简单的阈值法TSS浓度超出合理范围比如小于0或大于某个上限的像元直接置为NaN。虽然粗犷但很多L2产品本身就做过质量控制这一步只是兜底。经纬度和数据矩阵的维度关系也值得留意。GOCI Swath数据常见排列是列行或者行列读出来之后如果后续画图发现经纬网和浓度矩阵对不上很可能就是行列方向或者转置的问题。可以在单文件调试阶段就把size打出来对比确认不要等到批量出图才发现。4. 绘制TSS变化过程的实操4.1 单片快速预览与colormap选择TSS浓度图的核心是快速评估空间分布所以第一步先是单片渲染。我推荐pcolor加shading interp或者imagesc配合经纬度坐标。二者区别在于pcolor适合非均匀网格imagesc更轻快但坐标轴会整成网格索引需要额外处理。对于GOCI的标准产品经纬度网格通常是规则矩形用pcolor就行figure(Color, w); pcolor(lon, lat, tss_double); shading interp; axis image; colormap parula; caxis([0 200]); % 根据你的海区实际情况调整 colorbar; title([TSS , datestr(timeNum, yyyy-mm-dd HH:MM), (UTC)]);配色上需要注意TSS和叶绿素的呈现习惯完全不同。叶绿素常用深蓝到红光的多段渐变而TSS更适合用棕黄色系或土黄色调因为悬浮泥沙本身就给人一种“黄水”的直觉。Matlab自带的热色系colormap其实表现力一般我平时会用一个自定义的色带白色或浅蓝底代表低值黄色过渡棕红色代表高浓度。这样近岸高值区一眼就能看出来。4.2 多时相变化过程的拼图与动画表现GOCI一天内多个时相最简单的变化过程展示方式就是拼图。把同一天内按时间顺序排列的多幅图放进一个大figure里每幅图显示时间标注这样读者能直观看到潮周期内泥沙舌的推拉和扩散。拼图的核心是归一化。所有子图必须共用同一个颜色比例尺比如caxis固定为[0 200]否则不同时相浓度差异会被颜色拉伸掩盖。也不需要每张图都画colorbar只在第一张或最后一张保留即可避免视觉噪声。下面是一个典型的批量拼图骨架figure(Position, [50 50 1600 900]); for k 1:N subplot(3, 4, k); pcolor(lon, lat, tss_double); shading interp; axis image off; caxis([0 200]); colormap(map); title(datestr(meta(k).timeNum, HH:MM), FontSize, 10); end如果时相数超过十几个拼图排不下我会把图保存成序列帧再合成为GIF或视频。Matlab里用imwrite写GIF帧很常规但要注意每次都写缓存和延时参数我这里有一段常用的保存脚本figure; for k 1:N % ... 绘图代码 ... frame getframe(gcf); [A, map] rgb2ind(frame.cdata, 256); if k 1 imwrite(A, map, TSS_animation.gif, gif, Loopcount, inf, DelayTime, 0.5); else imwrite(A, map, TSS_animation.gif, gif, WriteMode, append, DelayTime, 0.5); end end这种方式做汇报演示非常好用领导或合作方不需要从头看十几张小图一个动图就懂了。4.3 特定点位与断面的时间序列曲线空间图解决“分布长什么样”的问题但很多人真正关心的是“某一个站位的TSS随时间怎么变”。这就需要从每个时相文件中提取对应像元的值再按时间排序画曲线。第一步是确定目标像元的行列索引。如果你能拿到目标站点的经纬度就用经纬度匹配targetLon 122.5; targetLat 37.0; [~, colIdx] min(abs(lon(1, :, :) - targetLon)); % 具体按矩阵维度调整 [~, rowIdx] min(abs(lat(:, 1) - targetLat));匹配时需要注意如果网格是列行结构索引顺序别搞反。还有一种办法是直接用鼠标在图上点选用ginput获取经纬度再回到批量循环里对应提取。第二步就是循环读取所有文件的TSS提取该像元的值tssSeries NaN(numel(meta), 1); for k 1:numel(meta) tss h5read(meta(k).path, tssPath); tssDouble double(tss) * scale offset; tssDouble(tss fillValue) NaN; tssSeries(k) tssDouble(rowIdx, colIdx); end figure; plot([meta.timeNum], tssSeries, -o); datetick(x, HH:MM); ylabel(TSS (g/m^3)); grid on;曲线出来之后配合前后潮汐数据或流向数据做交叉分析可以回答很多有意思的问题。比如涨潮期TSS是不是明显降低、落潮期是不是泥沙峰值出现、台风过境之后高浓度带有没有向远海扩散。5. 常见问题与排查技巧实录5.1 打不开文件或读取失败的几个诱因批量处理时最讨厌的就是循环到一半突然报错。我遇到过的第一类问题是h5read打不开文件报错信息五花八门但根源往往是文件本身不完整或者上传过程中损坏。这时可以用h5info单独测一下如果连h5info都失败基本就是文件坏了。第二类问题是路径太长导致Windows API拒绝访问。GOCI产品本身文件名就长再套上几层分类目录很容易超过路径长度限制。对策是把文件平铺在一个短路径目录下或者用“8.3短文件名”方式引用。第三类问题比较隐蔽文件句柄不释放。在循环里反复打开关闭是正常操作但如果你中途用h5create或h5write出了异常没能关闭句柄后续文件就可能一直被占用读取会报拒绝访问。调试期养成每读完一个文件就手动h5close.all的习惯能少掉一半莫名其妙的错误。5.2 数值异常与图像全黑全白图像全黑最常见的原因是没有把FillValue转成NaN直接拿整数最小值或哨兵值去画图整幅图的动态范围都被它拽走了。处理方式前面已经写了先根据属性把填充值替换成NaN再用isnan把无效区域盖成白色或透明背景。另一个常见问题是忘了乘比例因子。你如果发现图上的数值分布在“0到1”之间但实际量级应该在几十到几百那多半就是scale没乘。还有一种情况是数值分布正常但图面很花噪点很多这时要确认是否做了平滑。pcolor的shading interp有插值平滑作用但如果数据有异常的孤立小值平滑后周围会被拉出一片“墨渍”这时候要考虑对异常像元做中值滤波。5.3 坐标错位与显示变形经纬度矩阵和数据矩阵维度不一致是刚上手时最容易出的错。GOCI网格有时返回的是二维矩阵有时是三维含时间维的单帧如果你直接取lon(:,:,1)再和其他变量拼接维度就对不上。建议读取后立刻检查size(lon) size(tss_double)如果长宽刚好相反用permute或transpose处理不要用flipud硬翻因为经纬度和数据的方向关系在不同产品里并不一致。显示变形往往是不等经纬投影导致的。如果你直接plot经纬度数据而不设定axis equal近岸区域看起来会被拉伸。加上axis image或axis equal后能明显改善。需要的话还可以用m_map或自写墨卡托投影但基础分析阶段不需要过度投影直接坐标显示就够了。5.4 批量处理太慢怎么办大范围TSS数据动辄几千乘几千像元如果每天八个时相全读、全画Matlab很容易卡到怀疑人生。优化方向有几个按性价比排序第一只裁剪目标区域通过经纬度索引取出一个子矩阵再处理第二用parfor并行读取但要注意parfor里不变量的定义方式第三代码层面用matlabFunction或者预分配数组避免循环里动态增长。裁剪这招最粗暴也最高效。比如我只关心某海域一块区域先把经纬度范围的roi索引计算出来再在循环里只读那一块rowRange find(lat(:,1) minLat lat(:,1) maxLat); colRange find(lon(1,:) minLon lon(1,:) maxLon); tss_roi tssDouble(rowRange, colRange);注意裁剪必须在乘比例因子之前做这样处理的数据量直接小一个量级速度提升非常明显。6. 一些实操心得批量读取.he5这件事技术难度其实不算高真正吃时间的是跟产品结构较劲。我的经验是永远不要让脚本假设文件内部结构一成不变每一个新数据源到手第一件事就是跑一遍h5info。宁可多花五分钟看结构也好过后面输出一堆错误图再回头追查。另外TSS变化过程图的配色和比例尺一定要统一。很多人第一张图画得惊艳第二张图换了数据之后色标自动变了两张图放一起对比就完全失真。我在脚本里会把caxis写成一个全局变量所有时相都引用它从源头上避免这个问题。最后再分享一个小技巧GOCI产品里的经纬度矩阵和数据矩阵在对齐之后建议保存成一份.mat缓存文件。下次再读同一个时间段的数据直接load比逐文件h5read快得多。如果数据量比较大还可以只缓存ROI区域处理速度提升明显。这个习惯帮我度过了很多次deadline之前的批量重跑。