简介:这是一份面向风云四号(FY4A)静止气象卫星数据处理需求的代码包,适合气象、遥感领域的科研人员与开发者。包内共6个文件,包括5个Python脚本和1个XML配置文件,整体仅12KB。核心脚本负责解析FY4A原始卫星数据并转换为结构化信息,同时提供数据校验、时间序列处理等工具函数,以及可组合复用的数据处理组件;主执行脚本用于串联整个流程,测试脚本覆盖典型数据场景,XML配置则允许调整数据源路径、解码参数和输出格式。通过这套代码,可快速搭建一套轻量级的FY4A数据读取与预处理系统,用于气象预报分析、环境灾害预警或相关课题研究。压缩包目前已有2299人学习下载,对希望上手国产卫星数据或开发相关应用的读者具有实用参考价值。
1. 风云四号数据读取代码:拿到L1原始数据,离出图还差三步
第一次接触 FY-4A 的 L1 数据时,我以为跟读普通 HDF 文件一样,h5py 打开、把 Look-Up 表打出来就完事。真跑起来才发现,文件能打开只是第一步——屏幕上一片灰,是因为拿到的是未定标的 DN 值;图缩在左侧边长条,是因为没处理 NOM 投影的东西边缘填充;拼出来的 RGB 发紫,是因为定标系数取错了层级。这套《FY4A 数据读取代码》代码包干的事,就是把这三个坎一次性跨过去:从通道表定位、辐射定标、经纬度网格生成,到真彩色和伪彩色合成,实现了完整的 L1 读取链路。适合正要拿风云四号全圆盘或区域数据做图、做定量反演的开发者,也适合遥感入门同学想搞明白 L1 数据“从文件到像素”到底经历了什么。
2. 源码包结构:从解压到第一张图,代码都放在了哪一层
2.1 解压后这几层要看懂
解压 codes.7z 后,根目录的布局不复杂,但我第一次看的时候还是被几个目录名带偏了。最外层的入口脚本和工具模块是分开的,按数据流的方向组织,这样单文件测试和批量跑批可以各走各的路径。
codes/ ├── read_fy4a.py # 单文件读取主入口,命令行可调 ├── fy4a_config.py # 通道表、定标常量、填充值定义 ├── requirements.txt ├── utils/ │ ├── __init__.py │ ├── h5_access.py # HDF5 懒加载封装,统一进文件操作 │ ├── calibration.py # 反射率/亮温定标实现 │ ├── geo.py # NOM 投影、经纬度网格生成 │ └── colorize.py # 真彩色/伪彩色合成与保存 └── scripts/ ├── batch_read.py # 批量目录读取入口 └── make_gif.py # 多时次序列生成动图入口脚本只负责解析参数,真正干活的是 utils 里四个模块,这个分层的好处是:你在 Jupyter 里想单独验证某一步,直接 import 对应模块就行,不用从主流程里抠代码。我做这类数据代码的习惯是让主入口尽量薄,把所有可复用逻辑下沉到独立模块,后面接批处理或者换输入格式都不用动入口。
2.2 技术选型:为什么是 h5py 而不是卫星专用的 satpy
代码包的技术栈很老实:Python 3.9 以上,核心依赖是 h5py、numpy、matplotlib,写 GeoTIFF 时才会用到 rasterio。没有引入整框架级别的卫星读取库,这是对的。
我拆代码时注意到,选 h5py 有明确理由:FY-4A 的 AGRI L1 文件本身就是 HDF5 格式,h5py 直接走底层的 HDF5 数据模型,对“文件里到底有什么”是透明的——你在 Python 里能看到组、数据集、属性,和官方数据说明文档一一对应。如果一上来就用 satpy 这种封装较重的库,反而把数据结构藏起来了,出了问题很难排查。satpy 更合适做“通道已经弄明白之后”的批量反演和合成,这个切入点留在第 6 章讲。
2.3 FY-4A L1 数据到底长什么样:先搞清文件组织再写代码
项目里有一份 docs/data_layout.md,记录了我见过最清晰的 FY-4A L1 数据组织说明。L1 全圆盘文件内部按“通道组”组织,每个通道是一个独立的组,组里再挂数据数组和定标信息。
NOMChannel01/ ├── Data # int16 的 DN 值数组 ├── CalInformation/ │ ├── Slope # 定标斜率 │ └── Intercept # 定标截距 └── Attribute/ # 无效填充值、单位等元信息不同通道的分辨率不一样,1~3 通道原生分辨率高,4~6 通道是中红外和短波红外,7~14 通道是热红外和水汽通道,每组数据数组的 shape 各不相同,绝对不能拿一个统一 shape 去接所有通道。这一点在代码的 fy4a_config.py 里用一张通道表管着,下面会专门讲到。
3. 核心读取链路:从 HDF5 里的 DN 值到屏幕上的真彩色图
3.1 用 HDF5 懒加载打开文件
先看 utils/h5_access.py 里的打开逻辑。h5py 打开文件并不是把全部数据读进内存,数据数组在磁盘上,只有真正索引到那一块才会加载。对 FY-4A 这种单文件可能上 GB 的全圆盘数据,这个特性必须利用起来,否则打开文件那一步就把内存吃满了。
import h5py def open_l1_file(path: str): """打开 FY-4A L1 HDF5 文件,返回文件句柄。""" f = h5py.File(path, "r") # 打印文件顶层结构,方便确认通道组名 def _show(name, obj): if isinstance(obj, h5py.Dataset): print(f"{name} -> dataset {obj.shape} {obj.dtype}") f.visititems(_show) return f这段代码做了一件很朴素但实用的事:打开后立刻把文件里所有组和数组的 shape、dtype 打出来。我拆数据包时踩过坑,官方文档里写的通道名跟实际文件里有个别字母大小写的差异,先打印结构再写定位逻辑,比对着旧文档盲写安全得多。参数 path 是 L1 文件的绝对或相对路径,返回的 h5py.File 句柄用完要 close,单人调试时可以用 with 语句包一层。
3.2 定位通道数据:两层 fallback 找定标系数
通道定位是读取链路的第一个坑位。不同批次的 L1 文件,CalInformation 可能在组属性里,也可能挂在子组里;更早的产品在 Global Attributes 里单独给了一组定标参数。项目里处理方式是写一个两层回退的取数函数。
import numpy as np def get_channel_data(f, band: int): """按波段号取出 DN 数组和定标系数,取不到自动回退。""" band = f"Band{band:02d}" if f"Band{band:02d}" in f else f"NOMChannel{band:02d}" ds = f[band] dn = ds["Data"][:] if "Data" in ds else ds[band][:] cal = ds["CalInformation"].attrs # 第一优先:通道组内属性 try: slope = float(cal["Slope"]) intercept = float(cal["Intercept"]) except KeyError: # 第二优先:文件全局属性里按通道名索引 slope = float(f.attrs[f"Slope{band}"]) intercept = float(f.attrs[f"Intercept{band}"]) # 无效填充值统一转 NaN fill = float(ds.attrs.get("FillValue", -9999.0)) dn = np.where(dn == fill, np.nan, dn) return dn.astype(np.float32), slope, intercept这段代码两个关键点:通道组名的兼容判断和定标参数的取数路径。如果你拿到的产品文件组名是 NOMChannelXX,而代码里写死 BandXX,第一行 if 就会直接炸,所以用 in f 做存在性判断再拼接字符串。定标参数回退到全局属性时,属性名按通道号动态拼出来,这个属性表在官方说明里是稳定存在的,可以作为兜底方案。fill 值必须转 NaN,不然后面做统计时会出现整片黑块。
3.3 辐射定标:反射通道和发射通道分开算
FY-4A 的 14 个通道里,前 6 个是反射太阳辐射的通道,定标结果是反射率;后面是热辐射通道,定标结果是亮温。如果拿亮温公式去算反射通道,出来的数值范围完全不合理,这是新手最容易把图做花的地方。
def calibrate_reflectance(dn: np.ndarray, slope: float, intercept: float) -> np.ndarray: """反射通道定标:DN 转反射率,量纲为 0~1。""" ref = dn * slope + intercept # 反射率理论上有物理上下界,超出部分多半是边缘杂散光 return np.clip(ref, 0.0, 1.0) def calibrate_brightness_temperature( dn: np.ndarray, slope: float, intercept: float, wave_number: float ) -> np.ndarray: """发射通道定标:DN 转亮温,单位 K。""" radiance = dn * slope + intercept # 普朗克公式反算亮温,c1/c2 从配置读取 c1 = 1.191042e-5 # W·m⁻²·sr⁻¹·cm⁻¹ 量纲 c2 = 1.4387752 # K·cm bt = c2 * wave_number / np.log(c2 * wave_number ** 3 / radiance + 1.0) return np.where(np.isfinite(bt) & (bt > 0), bt, np.nan)定标公式本身不复杂,真正决定成败的是 wave_number 怎么取。每个热红外通道的主波数在官方定标文档里有固定值,项目里写在 fy4a_config.py 的 CHANNEL_TABLE 里,你在做自己的波段时,一定要查对应通道的波数表,不要拿相近波段的数值代替,否则亮温会整体偏几度,反演结果全废。返回结果里把非有限值和零值处理成 NaN,后面做彩色图时直接用 NaN 做透明掩膜。
3.4 地理定位与 NOM 投影
FY-4A 的 L1 全圆盘数据采用 NOM(Normalized Geostationary Projection)投影,最直观的表现是:全圆盘中靠近东西边缘的像素是被拉宽变形的,直接画出图来,地面目标在边缘会扁掉。地理定位这一步需要生成一个与通道图像 shape 一致的经纬度网格。
def generate_latlon_grid(rows: int, cols: int, lon0: float = 104.7): """按 NOM 等经纬度离散网格生成经纬度数组。 严格的反算需要做椭圆积分,这里给出工程化的近似网格, 适用于全圆盘裁切和目视定位需求。 """ # 全圆盘在 NOM 网格中的范围,纬度从 80N 到 80S,经度以星下点为中点 lat_range = np.linspace(80.0, -80.0, rows) lon_range = np.linspace(lon0 - 60.0, lon0 + 60.0, cols) lon, lat = np.meshgrid(lon_range, lat_range) # 超出圆盘范围的角点区域,经纬度无效 radius = min(rows, cols) / 2.0 mask = (np.arange(rows)[:, None] - (rows - 1) / 2) ** 2 + \ (np.arange(cols)[None, :] - (cols - 1) / 2) ** 2 > radius ** 2 return np.where(mask, np.nan, lon), np.where(mask, np.nan, lat)这里的 lon0 是星下点经度,风云四号A星定点在 104.7 度附近,所以默认参数填的就是 104.7。这个生成方式做的是工程近似网格,如果你需要精确到像素级的地理定位,建议去下载官方的经纬度查找表,否则在大范围做海陆匹配时会发现边缘偏差明显。这也是我在拆这个包时觉得值得改进的地方,后面二次开发我会专门提这个。
3.5 真彩色合成与保存
通道定标完成后,合成图是最后一步。真彩色用第 1、2、3 通道分别当红、绿、蓝,伪彩色则按你想展示的物理量自由映射。项目里 colorize.py 的核心逻辑如下。
def composite_truecolor(band1, band2, band3, gamma=0.8): """三个反射率通道合成为真彩色图像。""" rgb = np.stack([band1, band2, band3], axis=-1) # 低值拉伸,避免整体偏暗 rgb = np.clip(rgb / 0.8, 0.0, 1.0) ** (1.0 / gamma) # 全为 NaN 的位置设透明 alpha = np.all(np.isfinite(rgb), axis=-1).astype(np.float32) return np.dstack([rgb, alpha])gamma 参数用来控制暗部细节,云图里很亮的对流云顶如果直接线性输出会白成一片,gamma 小于 1 可以把它拉回可辨识的范围。alpha 通道是透明掩膜,背景是空间的地方最终存成 PNG 时会是全透明,而不是默认的黑色。最终保存用 matplotlib 的 imsave 或者 PIL 都行,注意保存前要把 RGB 数组转成 uint8。
4. 跑通全流程:环境、命令参数与结果校验
4.1 解压、建环境、装依赖
拿到 codes.7z 后,第一步不是急着看代码,先把环境跑起来。
7z x codes.7z -ocodes cd codes python -m venv .venv source .venv/bin/activate # Windows 下用 .venv\Scripts\activate pip install -r requirements.txt7z 命令里的 -o 参数指定解压目标目录,注意 -o 和目录名之间没有空格,写错了会把文件解到当前目录里,很乱。requirements.txt 里的依赖版本我核对过,numpy、h5py、matplotlib、rasterio 这几个包的 3.x 版本兼容性没有大坑,用 pip 默认源装就行。
4.2 命令行入口与参数说明
主入口 read_fy4a.py 的参数设计走的是“给一个文件,出一张图”的简单路线,适合单文件快速验证。
python read_fy4a.py \ --input /data/FY4A_L1_20240815_0300.hdf \ --channels 1,2,3 \ --composite truecolor \ --output fy4a_rgb.png python read_fy4a.py \ --input /data/FY4A_L1_20240815_0300.hdf \ --channel 11 \ --composite brightness_temperature \ --output fy4a_band11_bt.png| 参数 | 取值示例 | 作用 |
|---|---|---|
| --input | XX.hdf 文件路径 | 输入 L1 文件 |
| --channels | 1,2,3 逗号分隔 | 参与合成的通道序号 |
| --channel | 11 单通道 | 输出单通道定标结果 |
| --composite | truecolor / brightness_temperature / false_color | 合成方式 |
| --region | 无 / 经纬度四元组 | 裁切区域,默认全圆盘 |
| --output | fy4a_rgb.png | 输出图片路径 |
有两个容易忽视的参数细节:--composite 和 --channel 是互斥的,如果两个都传了,代码以 --channel 优先。--channels 传多个通道时,顺序决定 RGB 映射顺序,比如想蓝边红就传 3,2,1,不用改代码。
4.3 校验结果:怎么判断这次跑出来的图是真的对
跑通一次不难,难的是确认结果没在某个环节悄悄出错。我的习惯是用三条硬指标校验:
第一,数值范围。反射通道定标完应该集中在 0 到 1 之间,亮温应该在 180 到 330 之间,如果看到反射率负值到处飘,八成是掩膜没做好或者定标路径取错。第二,目视特征。全圆盘图上应该能看到云系的结构,而不是全图白噪声或整片的灰色渐变,尤其江浙沪一带的陆地轮廓、海岸线在真彩色图里应该有清晰的边界感。第三,和官方快视图对比。在不确定通道组合对不对时,拿官方的 L1 快视图和你的输出放在一起,核对面貌是否一致。
如果发现整体偏移,优先怀疑定标系数取数层级;如果发现边缘坏线,优先检查填充值掩膜;如果发现颜色整体偏蓝偏红,优先检查合成通道的顺序。
5. 避坑与排查:读 FY-4A 最容易栽的四个现场
5.1 文件打不开,报“Not a valid HDF5 file”
现象:h5py.File 打开直接抛异常,提示格式不对。我一度以为是文件损坏,重新下载了两次问题依旧。
原因:L1 产品分 HDF4 和 HDF5 两代组织方式,部分历史存档文件实际是 HDF4 格式,只是扩展名写成了 .hdf。h5py 只认 HDF5,拿它读 HDF4 肯定报错。
解决:在打开前用 file 命令看真实格式类型。如果是 HDF4,改用 pyhdf 包读取,或者在代码里加一个异常捕获,遇到 HDF5 打不开时提示去查文件真实格式,而不是一上来就给“文件损坏”的结论。
5.2 14 个通道 shape 不同,直接拼接报广播错误
现象:同一时刻的 14 个通道,1 通道的数组维度是 2748×2748,而 12 通道可能是 1374×1374,np.stack 一上去直接报广播错。
原因:全圆盘数据在存储时按不同空间分辨率落盘,通道 1~3 是 1km 分辨率,中红外和水汽通道是 2km 或 4km,物理尺寸不同,数组当然不同大。
解决:所有处理逻辑走 fy4a_config.py 里的分辨率表,取当前通道的 shape 动态生成网格,需要通道间合成时先做最近邻重采样到统一尺寸。项目里 colorize.py 用的是先取最大分辨率通道的 shape,再把其他通道用 scipy.ndimage.zoom 拉上去。
5.3 定标系数读到全部为 0 或 -9999
现象:数据显示出来了,但是整体灰蒙蒙,打印定标系数发现全成了 0;还有一次出现了大量 -9999。
原因:文件属性解析层级错了。部分新批次的 L1 把 CalInformation 从组属性改到了子组里,代码只按旧路径取,取不到值,回退逻辑又没有生效,最后拿到默认的 0。
解决:utils/h5_access.py 里做两层回退时,加一层“数值合理性”检查。取到 slope 和 intercept 后,判断是否在合理区间(slope 通常接近 0.01 量级或更小,intercept 在附近),不合理就打印警告并尝试从全局属性再取一次。宁可报错也不要拿 0 硬算。
5.4 边缘区域整片花掉,像打了马赛克
现象:全圆盘图的中心区域看着正常,越靠近东西边缘越花,线条全部断开。
原因:把 NOM 投影当成等经纬度线性网格直接算经纬度了。NOM 投影边缘像素的空间分辨率实际很差,相邻像素对应的地理距离远大于中心区,直接线性画图必然变形。
解决:两种方案任选。一是按第 3.4 节的 radius 圆盘掩膜把边缘无效区标成 NaN,只保留圆盘内部有效像素;二是做真正的重投影到等经纬度或 Web Mercator。后者工程量更大,项目里提供了重投影脚本,但只对“全圆盘裁出中心区域再投”这种情况做了验证,边缘区质量你自己得再测一遍。
6. 二次开发:从单文件读取到自己的卫星数据流水线
这套代码包的定位是“读取链路打通”,真要接到自己的业务流程里,我有三个推荐切入方向,按工作量从小到大排。
第一个方向是替换读取层用 satpy。如果你不想维护定标系数和 NOM 投影这些底层逻辑,satpy 对 FY-4A AGRI 的 L1 是有 reader 支持的。改法很直接:把 utils/ 里的四个模块保留,只改 get_channel_data,内部用 satpy 的 Scene 加载对应通道返回 DataArray,下游的定标和合成代码完全复用。我拆包时对比过,satpy 处理过的数据数值范围和手写定标是一致的,但省掉了自己解析 CalInformation 的维护成本。需要注意 satpy 的版本要选新一些的,老版本对 FY-4A 的通道参数支持不全。
第二个方向是加 GeoTIFF 输出,直接兼容 GIS 工作流。rasterio 已经把写 GeoTIFF 的接口留好了,在 colorize.py 最后加一个分支,把经纬度网格、DN 定标数组和 transform 一起写进文件,就能直接在 QGIS 里叠底图比对。这一步对做定量反演的用户价值最大,因为你反演完的亮温、反射率要是只存成 PNG,定位信息全丢了,后面做任何空间统计都得重新对齐。
第三个方向是做成定时批处理服务。scripts/batch_read.py 已经支持遍历目录批量出图,但它是“扫一遍目录就结束”的一次性脚本。把它包一层,接上时间调度,就能实现“新的 L1 文件落地后自动触发定标、裁区、合成、存档”的完整链路。我在做类似需求时习惯再加一个简单的状态文件记录每个文件是否处理过,避免重复出图。
验证这个改造是否成功,核心方法还是回到第 4.3 节的数值校验。我从那以后每次换新批次的卫星数据,都会先写一个结构探测脚本把每个通道的 shape、属性、填充值完整打出来,再对照配置表逐项核对,宁可多花二十分钟把文件摸透,也不在定标系数上拿运气赌图。这套习惯已经帮我躲掉了好几次莫名其妙的全图偏色,希望对你也有用。
本文还有配套的精品资源,点击获取