风云四号HDF5数据读取与处理：从标称投影到可视化全流程解析

简介面向风云四号卫星数据处理研究人员、气象遥感学习者以及相关应用开发者这份轻量级代码包提供了FY4A原始数据的完整读取与解析方案。压缩包共6个文件由5个Python脚本和1个XML配置组成整体仅12KB结构简洁便于快速定位核心模块。代码按职责划分核心模块将二进制卫星数据转换为结构化信息工具函数辅助数据校验、时间序列与预处理组件层封装图像生成和解码逻辑主程序统一调度流程测试脚本用于验证不同数据场景下的正确性XML配置则支持修改数据源路径、解码参数与输出格式。目前已有2299人学习参考尤其适合希望深入理解FY4A数据格式、自主构建处理流程的读者。借助这套代码可快速搭建从数据读取到分析的基础框架支撑气象研究、环境监测和灾害预警等领域的数据应用。1. 风云四号数据读取代码拿到HDF5文件的下一步怎么走不少同行第一次拿到风云四号FY-4A的L1、L2产品时第一反应是直接读文件结果要么是“不认识的HDF5结构”报错要么是读出来一堆无符号整型数完全不知道该往哪一步处理。风云四号的各类产品以HDF5为容器数据组织和常见的气象卫星格式不太一样如果不先了解它的层级关系写读取脚本时会频繁返工。这套阅读笔记所对应的代码包恰好把“打开文件、读数据集、按标准投影换算经纬度、把通道数据变成可视化影像”几个环节串成了可直接套用的脚本适合刚接触卫星数据的开发者和科研人员用最短的时间打通全流程。2. 前置理解FY-4A数据格式的层级与读取逻辑2.1 认识HDF5容器FY-4A产品的文件树结构风云四号产品文件通过HDF5Hierarchical Data Format version 5存储这类文件不像是平铺的二进制而是类似文件夹的树形结构。以AGRIAdvanced Geostationary Radiation Imager载荷的L1数据为例打开后最外层一般能看到//路径下挂载着Image、地理定位属性以及各通道的数据组。读取的第一步是搞清楚数据集的路径。代码包里提供了一个目录浏览函数作用就是让使用者快速列出文件内的所有组与数据集避免一次次盲猜路径import h5py def walk_hdf5_group(name, obj): if isinstance(obj, h5py.Dataset): print(fDATASET: {name} shape{obj.shape} dtype{obj.dtype}) else: print(fGROUP: {name}) def inspect_fy4a(file_path): with h5py.File(file_path, r) as f: f.visititems(walk_hdf5_group)这段代码利用visititems递归遍历文件内所有节点把所有数据集名称、形状和数据类型打印出来。逻辑上不需要对文件路径做任何假设适合拿到一份未知结构的产品时做侦察。参数方面如果文件较大打印输出会很长建议把结果重定向到txt文件例如在命令行追加 tree.txt。实际使用中会发现L1文件的通道编号一般以NOMChannelSunny、NOMChannelWV这类名称呈现而L2云产品通常是单独的数据组这种差异在刚开始接触时会觉得很乱跑一次目录树就会清晰很多。2.2 读取编码约定通道名称、填充值与数据维度不管后续要做什么处理读取之前必须确认三件事通道名称对应的物理含义、无效值使用了什么填充值FillValue、数据维度是行优先还是列优先。风云四号产品多数通道的维度是(行, 列)即第一维是南北方向第二维是东西方向这个约定与很多中等分辨率成像仪产品一致。常见的一个问题是不清楚有效数据范围。全套产品一般用无符号16位整型存储而无效或未观测区域的值会被设定为65335或65535这类特定数。如果在显示前不处理图像里就会出现大量白色噪点。先查看属性是稳妥做法with h5py.File(fname, r) as f: dset f[/NOMChannel01] attrs dict(dset.attrs) print(通道名称:, dset.name) print(填充值属性:, attrs.get(FillValue)) print(有效范围属性:, attrs.get(ValidRange))这里通过dset.attrs获取数据集的元数据能直接看到FillValue和有效范围。从实际经验看FY-4A L1数据有时也会把无效值信息写在文件级Attribute里并冠以CalibrationCoefficient之类的名字所以读数据集前先打印attrs是省时间的关键动作。做完这一步后续的掩膜和定标就不容易翻车。3. 核心读取流程从标称投影到经纬度反算3.1 地理定位用NORTH_EAST坐标计算逐像素经纬度风云四号L1全圆盘数据采用标称投影文件里没有直接提供经度和纬度矩阵而是通过一个“地球固定网格”换算过来。换算需要用到文件属性里的中心经度、比例因子等参数。实际项目中最常用的做法是采用解析几何法先计算扫描角再转换为经纬度。import numpy as np def compute_latlon(image_shape, center_lon104.7): rows, cols image_shape # 全圆盘标称网格的步长单位为弧度 step (0.5 / 3600.0) * np.pi / 180.0 # 行列索引转成相对中心的角度 col_angle (np.arange(cols) - (cols - 1) / 2.0) * step row_angle (np.arange(rows) - (rows - 1) / 2.0) * step # 这里采用近似共形投影反解 lon center_lon np.degrees(col_angle) lat np.degrees(row_angle) lon2d, lat2d np.meshgrid(lon, lat) return lon2d, lat2d这段代码实现了最简版本的地理定位换算。关键参数是step全圆盘图像1像素对应0.5角秒转弧度后用于计算每个像素相对星下点的偏移。如果对精度要求更高则需要引入地球椭球半径与卫星高度进行迭代反算但大多数可视化场景下用这个近似版本已经足够。使用的时候需要注意meshgrid的维度顺序要和图像数据的行、列对应否则画出来的等值线会整体错位。这个近似方案在早期调试阶段特别有用能快速验证数据是否读取正确。3.2 通道数据预读DN值转反射率与亮温对于可见光通道L1产品提供的是计数值DN值需要结合定标系数转成反射率。L1文件里一般会附带各通道的定标斜率Scale和截距Offset读取后直接用广播方式计算即可def dn_to_reflectance(dn, valid_min0, valid_max4095): # 典型FY-4A可见光通道定标参数会写在文件attrs里 scale 0.001 # 示例值实际以文件为准 offset 0.0 reflectance dn * scale offset # 无效值处理超出有效范围的置为NaN reflectance np.where((dn valid_min) (dn valid_max), reflectance, np.nan) return reflectance红外通道则通常使用普朗克公式的反函数把辐射值转为亮温工程上常用查表法。查表法的主要思路是针对通道建立DN到亮温的一一映射这种方式的好处是计算速度快实现时只需构造一个数组而不用逐个点调用数学函数。参数层面的两个重点是有效范围必须依据文件属性来确定以及确认定标系数是文件给定。之前有同行因为套用可见光通道的定标公式读取红外数据结果亮温起点差了十几K费了很长时间才定位到问题。真实使用中建议保留原始DN值在内存中不要就地覆盖因为后续做云检测或质量控制时有些算法需要原始计数值重新计算。4. 分类型读取L1辐亮度与L2云产品的不同处理4.1 L1辐亮度产品通道合并与黑白图像导出L1产品是所有应用的基础最常见的操作是从多通道数据中选取若干通道合成RGB图像。FY-4A的AGRI载荷的可见光通道在晴空条件下特征清晰直接取三个不同谱段的通道作为R、G、B分量进行简单拉伸后可以得到类似真彩色的效果。注意风云四号并没有和卫星平台其他载荷完全一样的真彩色三通道命名需要根据各通道的中心波长自行映射组合。from PIL import Image def scale_to_uint8(data): data np.nan_to_num(data, nan0.0) data (data - np.nanmin(data)) / (np.nanmax(data) - np.nanmin(data)) return (data * 255).astype(np.uint8) def save_rgb_from_channels(ch1, ch2, ch3, output_namergb.png): r scale_to_uint8(ch1) g scale_to_uint8(ch2) b scale_to_uint8(ch3) rgb np.stack([r, g, b], axis-1) Image.fromarray(rgb).save(output_name)这段代码演示了合成输出。核心是把不同通道的DN值归一化到0到255的整数范围再用np.stack把三个二维数组合并为三维矩阵。实践中需要注意三个通道各自的动态范围可能差别很大直接把各自数值做全局最小最大归一化会导致某一通道在合成图上被压暗。此时需要对每个通道做百分位截断比如保留2%到98%之间的数据。代码包里如果只做了基础归一化建议自行加上分位数裁剪视觉上会更接近传统效果。4.2 L2云产品读取分类编码与有效值掩膜云检测、云相态这类L2产品与L1不同不再使用辐亮度作为数值输出而是用整数编码代表类别。例如云掩膜产品里0一般代表晴空1代表有云2代表不确定或云影边缘。读取时要格外注意类别码与填充值的区分def read_cloud_mask(fname, dataset_path/Cloud_Mask): with h5py.File(fname, r) as f: raw f[dataset_path][:] fill f[dataset_path].attrs.get(FillValue, 255) # 掩膜无效值 mask np.where(raw fill, 1, 0) # 返回原始数据和有效标志 return raw, mask这里的raw是原始整数编码mask是二值有效数据标志。后续做统计时一般只对mask 0的像素做类别频率分析这样不会把填充区域当成云。对于目录下多个文件需要做时间序列统计的场景读取端最值得优化的就是这类掩膜操作尽量用numpy的向量化比较替代循环。5. 避坑排查风云四号数据读取高频错误排查手册5.1 文件能打开但打印速度极慢现象使用visititems扫描文件结构时脚本卡住迟迟不结束。 原因L1全圆盘文件的数据集数量多部分数据集通过压缩存储遍历时触发解压。 解决不要对全文件遍历直接在已知路径下读取目标数据集。如果不知道自己需要的东西在哪个路径先通过h5py.File(fid, r)打开后用f.keys()查看顶层组再逐层进入。5.2 图像输出全黑或者全是白点现象导出的图像在图像查看器中是一整片白色或者整片黑色。 原因大概率是填充值没有正确处理。未观测区域的DN值较大直接参与缩放后被拉伸为白色或者选择的最小值包含了无效值导致有效像素被压缩到无法分辨的区间。 解决在归一化前先把无效值设为NaN然后再做分位数截取。截取最好在2%和98%分位处做必要时用直方图确认数据范围是否合理。5.3 经纬度定位和地图叠加时整体偏移现象把计算得到的经纬度矩阵叠加到海岸线上时发现整幅图偏移了十几公里。 原因中心经度参数未从文件属性读取。常见做法是硬编码104.7度但部分区域产品例如区域拼图的中心经度可能不同。 解决检查文件属性中是否包含Center_Longitude或者geostationary projection origin字段用属性值替换代码里的硬编码常量。通过属性动态获取会明显减少返工。5.4 红外亮温出现明显跳变的条带现象亮温图上出现横向或纵向的细密条带且不随缩放比例消失。 原因L1数据中部分扫描带存在坏元或噪声书面叫法是噪声条带。很多L1产品未对这种条带做完全修复。 解决对通道数据做中值滤波或者在合成前把该条带对应的行、列直接掩膜。需要根据图像上条带的具体位置和宽度写掩膜数组不建议用全局平滑否则会损失细节。5.5 16位整型数据转浮点时值域被截断现象读取后的数值最大只到65535但处理时发现有些负值或小数丢失。 原因HDF5中存储的数据类型是uint16如果直接用astype(np.float16)或默认的Int类型转换会丢失NAN的表示能力。 解决先转为np.float32再做定标计算。在需要存储中间结果时用np.float32这样既能保存小数也能通过设置NaN保持无效值语义。5.6 多文件读取时文件句柄泄漏现象连续读取几百个文件后内存占用持续上升处理速度越来越慢。 原因由于打开了h5py.File但未关闭导致文件描述符耗尽。 解决统一使用with h5py.File(...) as f:上下文管理器同时不要在循环内部保存大数组引用。读取完立即对数组做提炼运算临时变量用del后再下一次循环内存占比才稳得住。6. 进阶用法把读取脚本变成批量处理工具6.1 遍历目录中的多时间文件并提取特征单个文件能正常读取后大部分人的下一步是会获取某个时间范围内的连续观测数据用来做时序分析或制作动画。批量处理的框架并不复杂关键在于把文件路径的解析与数据读取解耦import glob from pathlib import Path def process_single_file(fname): with h5py.File(fname, r) as f: ch1 f[/NOMChannel01][:] ch2 f[/NOMChannel02][:] # 做一个通道差的指数作为简易云指数示例 index (ch1.astype(float) - ch2.astype(float)) / (ch1.astype(float) ch2.astype(float) 1e-6) return np.nanmean(index) def batch_process(root_dir): files sorted(glob.glob(str(root_dir / *.HDF))) results [] for fname in files: try: results.append((Path(fname).name, process_single_file(fname))) except Exception as e: print(f跳过文件 {fname}: {e}) return results这段代码中batch_process先使用glob模块匹配所有HDF文件排序后逐个处理。单文件异常通过try-except捕获不会因为个别坏文件中断整个批处理。对气象数据的批处理任务这一点非常实用文件多、存储介质差异大总会有某个文件损坏或拷贝不完整加一层异常隔离能省掉人工查找的麻烦。索引计算中用1e-6避免分母出现零值属于常规稳定化处理。6.2 把定位结果与通道数据一起写成GeoTIFF很多应用场景要求把处理结果交给GIS工具加载最直接的做法是写GeoTIFF。实现时不需要自己实现GeoTIFF编码器直接借助已有的osgeo.gdal库完成pip install GDALfrom osgeo import gdal, osr def write_geotiff(output_path, data, lon, lat): rows, cols data.shape driver gdal.GetDriverByName(GTiff) out_ds driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) out_ds.GetRasterBand(1).WriteArray(data) # 定义WGS84投影 srs osr.SpatialReference() srs.ImportFromEPSG(4326) out_ds.SetProjection(srs.ExportToWkt()) out_ds.SetGeoTransform([lon.min(), (lon.max()-lon.min())/cols, 0, lat.max(), 0, (lat.min()-lat.max())/rows]) out_ds.FlushCache()写入GeoTIFF的要点是GeoTransform六参数第一项是左上角经度第二项是像素宽度第四项是左上角纬度第六项为负值表示行向下。如果直接把lat矩阵的最小值和最大值放进去要保证lat的计算方向与数组排列一致否则输出的图像会上下颠倒。写完后用GIS软件叠加检查一次如果发现翻转把第六项和第四项对调即可。批量工具稳定后每天的数据处理基本就是一条命令的事。运行脚本前先对一两个文件做单文件测试测试时把中间变量存成NumPy的.npz格式既方便断言数值范围又便于后续直接复用处理中间结果。那以后我每换一批数据都会强制走一遍“目录树检查—填充值确认—单文件可视化—批量运行”这四步虽然多花几分钟但基本杜绝了临时改脚本的尴尬。希望这份从格式认知到批量处理的笔记能帮你在风云四号数据处理上少走一些弯路。本文还有配套的精品资源点击获取