MATLAB批量读取GOCI .he5文件提取TSS浓度与可视化详解

上周帮一位同行排查了一个很典型的痛点从GOCI水色数据的.he5文件里提取TSS浓度还要把连续几十个时相的悬浮物变化过程画出来。他用MATLAB一个个文件手动读读得人都麻了——这不怪他.he5这种HDF5封装格式对不熟悉的人确实不友好。但一旦摸清了读取路径和批量循环的思路整套处理流程其实可以做得非常顺滑。这篇文章我就围绕这件事展开怎么用MATLAB批量读取GOCI的.he5文件怎么把里面的TSS数据提取、还原并整理成时间序列最后又如何把变化过程用静态图、时序曲线和动态影像展示出来。适合刚接触GOCI数据处理、或者对HDF5格式不太熟又被批量这两个字卡住的人参考。我会把每一步为什么这么做的逻辑也讲清楚而不是只丢给你一段能跑的代码。1. GOCI水色数据与.he5文件先搞明白你要拆的是哪种数据1.1 GOCI卫星和TSS反演的基本逻辑GOCI是静止轨道海洋水色成像仪有8个可见光到近红外波段空间分辨率大概500米最大的特点是时间分辨率高——一天之内能对同一片海域进行多次观测几乎做到逐小时成像。这个特性对TSS总悬浮物浓度监测非常关键因为悬浮物受潮汐、风浪、河流输沙影响很大可能一两个小时就变一个样。极轨卫星一天只能过境一两次根本抓不住这种快速变化而GOCI能提供连续的高频观测序列这也是大家愿意花力气处理它的根本原因。TSS是水质评价里的核心参数通常用单位体积水体里悬浮颗粒物的质量来表示。在海洋水色遥感中它并不是直接测量出来的而是通过水体反射率信号反演得到的。水体里悬浮物越多后向散射越强卫星看到的反射率就会在某些波段明显升高。所以GOCI的L2级产品除了反射率数据外往往会附带一系列反演产品TSS就是其中之一。1.2 .he5文件的内部结构解剖.he5这个后缀名其实指的就是HDF5格式不过是官方在分发数据时约定俗成用.he5来标记。很多人的第一反应是找不到对应工具实际上MATLAB的h5系列函数就能直接处理完全不依赖第三方插件。第一次拿到文件别急着写代码先花一分钟用h5disp看看里面到底有什么h5disp(COMS_GOCI_L2_GA_20200101_000000.he5);你会看到一组嵌套的组和数据集典型的路径长这样/HDFEOS/GRIDS/GOCI_Grid/Data Fields/TSS最前面的是HDFEOS这一层固定结构中间是网格名称真正存数据的地方在Data Fields组下面。除了TSS这里通常还会有Rrs各波段、Lw各波段、QualityFlag等数据集。文件名里的GA一般表示这是地理投影产品坐标已经规则化可以直接按行列索引。我建议第一次拿到文件时把h5info的结果也打出来看一眼尤其是各个数据集的Attributes。因为很多关键信息比如无效值、缩放系数、偏移量都藏在属性里光看数值本身是看不懂的。1.3 TSS数据到底是怎么来的GOCI的官方L2产品在不同版本里TSS字段不一定都有。有的版本直接提供TSS产品有的版本只给遥感反射率Rrs需要自己反演。如果你拿到的文件里有TSS字段读取就直接用如果没有就得用Rrs波段自己算。常见的悬浮物反演思路是利用可见光红波段和近红外波段的反射率组合因为这两个波段对悬浮物浓度变化最敏感。比较常见的形式是波段比值TSS A * (Rrs(red) / Rrs(NIR)) B也有用对数形式的TSS 10^(a b * log10(Rrs(red) Rrs(NIR)))但这里的A、B或a、b绝对不是普适常数不同水体的光学属性差异很大系数必须用现场实测水样做区域标定。如果你只是为了看相对变化趋势可以用文献里公开发表过的系数做个近似如果要报定量结果必须重新拟合。2. 环境准备MATLAB读取he5文件的三种思路与我的选型逻辑2.1 主力方案h5disp h5read原生读取如果你用的MATLAB版本是R2016b之后推荐直接用h5系列函数核心就四个h5disp打印文件结构h5info获取结构信息并存入变量h5read读取数据集数据h5readatt读取属性值一个最基础的读取流程长这样filename COMS_GOCI_L2_GA_20200101_000000.he5; tss_path /HDFEOS/GRIDS/GOCI_Grid/Data Fields/TSS; tss_raw h5read(filename, tss_path); scale h5readatt(filename, tss_path, scale_factor); offset h5readatt(filename, tss_path, add_offset); tss double(tss_raw) * scale offset;需要注意老版本MATLAB里有hdf5read系列函数但从某个版本开始推荐用h5系列。新代码里直接用h5read就好别再依赖hdf5read因为官方对后者的维护已经不如前者积极。2.2 辅助方案用GDAL转一手数据有些场景下官方提供的.he5文件在MATLAB中读取会遇到版本兼容问题比如HDF5的某些复杂数据类型在MATLAB里支持得不好。这时候我一般会借助GDAL把数据先转成GeoTIFF或NetCDF再读入MATLAB。GDAL对HDF5的支持非常成熟尤其是各种地理参考信息都能正确识别。gdalinfo COMS_GOCI_L2_GA_20200101_000000.he5 gdal_translate -of GTiff HDF5:COMS_GOCI_L2_GA_20200101_000000.he5://HDFEOS/GRIDS/GOCI_Grid/Data Fields/TSS tss_20200101_000000.tif转换的好处是后续可以用readgeoraster直接读而且GeoTIFF自带地理坐标画图时叠加海岸线、经纬度网格都方便很多。坏处是多了一道工序批量处理时磁盘占用会增加不少。所以我通常只在单景调试或MATLAB读取报错时才用GDAL日常批量处理还是走纯MATLAB路线。2.3 选型逻辑为什么原生读取优先我用过的方案对比下来结论很明确优先用原生h5读取只有在文件格式异常时才考虑GDAL转换。原因有三个。第一GOCI数据的.he5文件内部就是标准的HDF5简单数组不涉及复杂嵌套类型MATLAB的h5read处理这类数据非常稳定完全没必要中间转一道。第二批量处理时延展性更好代码逻辑统一出现问题时也容易定位。第三h5read支持直接读取某个局部区域如果你只关心某个子区域可以先读坐标数组再切片避免整景数据全部载入内存。从实际工程角度看能少一道转换就少一道永远是省时省事的原则。3. 批量读取he5并整理TSS的完整操作流程3.1 拿到文件清单并按时间排序批量处理的第一步不是写循环而是先把文件列表整理干净。GOCI的L2文件名通常自带观测时间格式一般像这样COMS_GOCI_L2_GA_20200101_000000.he5 COMS_GOCI_L2_GA_20200101_010000.he5中间的日期和时间串可以直接解析。我习惯直接用dir列出所有he5文件再用正则表达式或简单索引把时间信息提取出来转成datetime数组然后排序。这一步看似简单但如果你跳过它后面画时间轴的时候绝对会乱。files dir(*.he5); fileNames {files.name}; dates NaT(numel(fileNames), 1); for i 1:numel(fileNames) tk regexp(fileNames{i}, (\d{8})_(\d{6}), tokens); dates(i) datetime([tk{1}{1} tk{1}{2}], InputFormat, yyyyMMddHHmmss); end [dates, idx] sort(dates); fileNames fileNames(idx);这里的datetime是绝对时间后面画横轴、标时间都用它比把时间当字符串处理方便得多。3.2 循环读取核心数据集有了文件列表循环读取的主体就简单了。但有几个细节我要强调一下第一空间参考信息经纬度数组或行列坐标每景都一样应该在循环外只读一次不要每景都重复读第二TSS数值最好用single存不要用double否则数据量直接翻倍第三用try-catch包一下读取过程个别文件损坏时不要中断整个循环。nFiles numel(fileNames); info h5info(fileNames{1}); % 假设经纬度数组在经纬度数据集里 lon h5read(fileNames{1}, /HDFEOS/GRIDS/GOCI_Grid/Data Fields/lon); lat h5read(fileNames{1}, /HDFEOS/GRIDS/GOCI_Grid/Data Fields/lat); [rows, cols] size(lon); tssCube zeros(rows, cols, nFiles, single); for i 1:nFiles try tssPath /HDFEOS/GRIDS/GOCI_Grid/Data Fields/TSS; tssRaw h5read(fileNames{i}, tssPath); s h5readatt(fileNames{i}, tssPath, scale_factor); o h5readatt(fileNames{i}, tssPath, add_offset); tssCube(:, :, i) single(tssRaw) * s o; catch ME warning(读取失败%s原因%s, fileNames{i}, ME.message); end end如果你的数据是原始TSS浮点产品没有scale_factor和add_offset属性那直接cast成single就行。关键是先看属性再决定要不要缩放千万不能想当然。3.3 还原真实TSS值与质量掩膜从.he5里读出来的原始值绝大多数时候不是真实的物理浓度而是一个编码后的整数值。GOCI的L2产品为了降低存储空间通常把浮点TSS乘上一个缩放系数再存成整数。如果你直接把整数当成真实浓度去画图数值差好几个数量级画出来的图完全不能用。同时还要处理无效值。卫星数据的无效值一般有两种一种是填了指定填充值比如-9999另一种是质量标记不通过的区域比如云覆盖、陆地掩膜。正常做法是读取时把这些位置的TSS值统一置为NaN这样后面的统计和绘图都会自动忽略这些位置。fillValue h5readatt(filename, tssPath, _FillValue); tssReal(tssRaw fillValue) NaN; % 如果还有质量标记一般大于等于某个阈值才视为有效 qf h5read(filename, /HDFEOS/GRIDS/GOCI_Grid/Data Fields/QualityFlag); tssReal(qf 64) NaN; % 这里的64只是示例实际阈值看产品说明这里要特别说一下阈值只是示例这个问题。质量标记的有效阈值因产品版本和机构而异不要拿别人的代码里的数字直接套用。最靠谱的做法是读取后把QualityFlag的分布直方图画出来先弄清楚0、1、2这些标记都代表什么再定阈值。3.4 构建统一的时间序列矩阵批量读取的最终产物我建议整理成一个三维矩阵tssCube维度分别是行、列、时相。这个结构的好处是后面无论想做空间分布、时间剖面还是逐点曲线都可以很方便地切片。% 某一点的时间序列 pointTss squeeze(tssCube(r0, c0, :)); % 某一时相的空间分布 snapshot tssCube(:, :, i); % 全区域逐时相均值 seriesMean squeeze(mean(tssCube, [1 2], omitnan));这一步做扎实之后后续所有分析都变成对tssCube的切片操作代码会简洁很多。我一般还会顺手把tssCube、lon、lat、dates一起存成.mat文件后面画图直接load不用重新读he5。如果磁盘空间不是问题也可以每个时相单独保存一个tif画图时按需读取。4. TSS变化过程可视化静态图、时序曲线与动态影像4.1 单时相TSS分布图的基本画法拿到数据后的第一个常规操作是把某个时相的TSS空间分布图画出来看看整体形态是否合理。GOCI的GA产品已经做过地理投影经纬度数组可以直接用来做坐标轴标签。figure; pcolor(lon, lat, tssCube(:, :, 1)); shading flat; colormap(parula); colorbar; caxis([0 50]); % 具体范围根据你研究区的实际情况调整 xlabel(经度); ylabel(纬度); title(TSS分布 - 第一时相);用pcolor而不是imagesc是因为pcolor能直接处理经纬度坐标不会出现行列数与经纬度不匹配的问题。shading flat是为了去掉网格线让图面干净。颜色轴的范围非常重要先看数据整体分布再定别一开始就凭感觉写死不然很容易把细节压没了。4.2 多时相变化过程固定点位的时序曲线TSS变化过程最直观的展示方式之一是选几个感兴趣的位置比如河流入海口、近岸观测站、某个羽流扩散带画出它们随时间的浓度变化曲线。这个操作的代码很简单关键是点位怎么选。% 假设r0,c0是第一个点位r1,c1是第二个点位 curve1 squeeze(tssCube(r0, c0, :)); curve2 squeeze(tssCube(r1, c1, :)); figure; plot(dates, curve1, o-, LineWidth, 1.5); hold on; plot(dates, curve2, s-, LineWidth, 1.5); legend(站点1, 站点2); xlabel(时间); ylabel(TSS浓度);选点位的时候建议先在单时相图上用数据游标看一下具体坐标对应的地理位置避免点选到陆地或始终有云覆盖的区域。潮汐作用明显的海域TSS曲线通常会呈现周期性涨落这时候时间轴千万别弄错细节我在后面第五节再展开。4.3 动态影像展示把变化过程做成动画比静态序列图更有冲击力的是把连续时相的TSS分布做成动画。一天之内你可以看到悬浮物羽流如何随潮汐进退这种动态表达在研究报告里说服力很强。我一般优先输出GIF因为通用性好PPT和文档里都能直接插入。用imwrite循环写帧就行figure; clim [minTSS, maxTSS]; for i 1:nFiles pcolor(lon, lat, tssCube(:, :, i)); shading flat; colormap(parula); caxis(clim); title(datestr(dates(i), yyyy-mm-dd HH:MM), Interpreter, none); % 写入GIF frame getframe(gcf); [im, map] rgb2ind(frame.cdata, 256); if i 1 imwrite(im, map, tss_change.gif, LoopCount, Inf, DelayTime, 0.5); else imwrite(im, map, tss_change.gif, WriteMode, append, DelayTime, 0.5); end end这里最关键的一行是caxis(clim)。如果每帧都让MATLAB自动调整颜色轴浓度高的时相会把整体色标拉高浓度低的时相又会变暗前后帧之间颜色根本无法对比。只有固定了颜色轴范围动画才能真实反映浓度相对变化。如果想输出MP4就把getframe的帧输入VideoWriter原理一样。4.4 可视化中容易犯的几个错误第一线性色标下低浓度区细节全丢。很多近岸水体TSS的动态变化其实发生在低浓度端如果你用线性色标画低值区全是蓝色一片。这种情况下我建议做一个log变换显示[~, ~, imInd] histeq(...)太复杂了直接用log10(tss)当色标输入就好颜色条上标注原始值即可。第二只展示最大值时相而不展示变化幅度。变化过程比绝对值更有信息量。我通常会额外画一张每个像素的TSS标准差图或最大差值图一眼就能看出哪里变化最剧烈这个图往往比单时相分布图更能说明问题。第三文件多的时候别所有图都存成一个巨大的figure。一次性plot几十条曲线没有任何可读性我更倾向于分组展示或者用子图把有代表性的几个时相并列。5. 实操中踩过的四个坑与对应解法5.1 坑一数据字段名和文件结构随版本变化第一次跑通之后你以为以后都简单了结果换一批数据突然报错提示找不到TSS数据集。这就是GOCI产品的结构在不同版本之间不完全一致造成的。有的版本把TSS放在Data Fields下有的版本可能嵌套更深有的字段名加了后缀。解决这个问题我一般写一个小探测函数自动遍历整个HDF5结构把所有数据集路径全部打印出来function listDatasets(info, prefix) if isfield(info, Datasets) for i 1:numel(info.Datasets) fprintf(%s/%s\n, prefix, info.Datasets(i).Name); end end if isfield(info, Groups) for i 1:numel(info.Groups) listDatasets(info.Groups(i), [prefix / info.Groups(i).Name]); end end end然后对任意一个.he5文件调用这个函数5秒内就能搞清内部结构。这个思路比死记住某个固定路径可靠得多强烈推荐你把这个函数存下来以后遇到任何HDF5格式的数据都能用。5.2 坑二无效值、负值、离群值没处理干净GOCI的TSS反演产品有个老毛病近岸高浊度区域偶尔会出现负值云边缘和水体很清的区域偶尔会出现离谱的大值。如果你直接画图色标会被这些异常值拉变形像样图里的细节全被淹没。我的处理习惯是读取后先做一次物理范围过滤然后统一置NaN。每个区域的水体条件不同合理范围也不同但你可以先用百分位数快速看一下数据的统计分布tmp tssReal(:); valid tmp 0 isfinite(tmp); prc prctile(tmp(valid), [1, 99]); tssReal(tmp 0 | tmp 10*prc(2)) NaN;这里的10倍上限只是一个粗略的离群值截断具体你得结合研究区的实际情况。重点是在处理流程中加入这个过滤步骤能避免后续画图和统计被极个别像素带偏。5.3 坑三内存占用和批量处理速度GOCI单景数据如果按2780乘以2780计算每个时相用single存大概是30多MB几十景加起来就是1GB上下的三维矩阵。看起来不大但如果你用double内存直接翻倍加上中间变量稍微不留神就能把普通办公电脑的内存吃满。实际处理建议就三条。第一三维矩阵一律用single不要贪方便用double。第二经纬度、掩膜这类空间参考只读一次固化到变量里循环里绝不重复读取。第三中间结果及时落盘比如把tssCube保存成.mat或逐时相tif后面画图时用多少读多少而不是每次从头跑一遍批量读取。如果你有并行计算工具箱循环体里用parfor确实能省时间但要注意GOCI文件分布在网络盘的话并发读取反而会更慢瓶颈在磁盘I/O而不在CPU。5.4 坑四时间标准和时区搞错这是最隐蔽也最致命的坑。GOCI文件名里的时间通常是UTC而你所处的研究区域用的可能是UTC9之类的地方时。如果你直接把文件名时间当成地方时来画时间轴TSS的变化峰值会整体偏移好几个小时恰好和潮汐相位对不上你还会以为是反演算法出了问题。我的做法是解析文件名时间之后统一在代码里显式标注时区。比如用datetime的TimeZone属性统一转成UTC或地方时后续所有分析都基于同一个时间标准dates datetime([...], InputFormat, yyyyMMddHHmmss, TimeZone, UTC); datesLocal dates; datesLocal.TimeZone Asia/Shanghai; % 按需设置成你研究区的时区这样至少在图上标注的是统一的时间不会出现同一批数据里一部分按UTC画、一部分按本地时间画的错误。另一个和时间相关的细节是GOCI一天的观测次数有限相邻两景之间可能隔一两个小时作时序曲线时最好把缺失时相作为断点处理别把不连续的时间点硬连成连续的线。最后想分享的一点经验整套流程跑下来我个人体会最深的不是某个函数怎么用而是先探测结构、再批量处理这个顺序千万不能省。我第一次写批量循环的时候就是拿到文件直接套路径读取结果换了数据版本后跑了一整夜早上一看80%的文件都读出来是NaN浪费了时间不说还差点把原始数据覆盖了。后来养成的习惯是无论数据换了哪个版本先花五分钟用h5disp或那个递归函数看一眼结构确认数据字段、属性、缩放系数都存在再启动批量流程。另外一个小技巧是批量处理前先挑一景数据从头到尾跑一遍包括读取、还原、画图、保存确认每一步的输出都符合预期后再动用循环。看起来多花了一点时间实际上帮你省掉的返工时间绝对值得。GOCI这种高频水色数据在近岸水质监测里潜力很大把读取这关打通之后你会发现后面不管是做潮周期分析、羽流扩散追踪还是悬浮物输运通量估算都只是从tssCube里切数据的事。