企业官网建设流程全解析

简介：面向风云四号卫星数据处理研究人员、气象遥感学习者以及相关应用开发者，这份轻量级代码包提供了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产品的文件树结构

风云四号产品文件通过HDF5（Hierarchical Data Format version 5）存储，这类文件不像是平铺的二进制，而是类似文件夹的树形结构。以AGRI（Advanced Geostationary Radiation Imager）载荷的L1数据为例，打开后最外层一般能看到//路径下挂载着Image、地理定位属性以及各通道的数据组。

读取的第一步是搞清楚数据集的路径。代码包里提供了一个目录浏览函数，作用就是让使用者快速列出文件内的所有组与数据集，避免一次次盲猜路径：

import h5py def walk_hdf5_group(name, obj): if isinstance(obj, h5py.Dataset): print(f"DATASET: {name} shape={obj.shape} dtype={obj.dtype}") else: print(f"GROUP: {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_lon=104.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_min=0, valid_max=4095): # 典型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, nan=0.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_name="rgb.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 GDAL

from 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格式，既方便断言数值范围，又便于后续直接复用处理中间结果。那以后我每换一批数据都会强制走一遍“目录树检查—填充值确认—单文件可视化—批量运行”这四步，虽然多花几分钟，但基本杜绝了临时改脚本的尴尬。希望这份从格式认知到批量处理的笔记，能帮你在风云四号数据处理上少走一些弯路。

本文还有配套的精品资源，点击获取