葵花8/9卫星数据处理实战:从HSD解析到云图动画生成
第一次被安排把葵花8的云图动画拉出来的时候我以为下载几张图拼一拼就能交差。等真正打开下载目录才发现里面躺着的全是.DAT结尾的二进制文件别说云图连图片头都找不到。那段时间我一边翻文档一边调程序从数据下载、HSD 格式解析、辐射定标、真彩色合成到区域裁剪和动画输出把葵花8/9 的完整处理链路老老实实走了一遍。后来这套流程成了我处理静止轨道卫星数据的基本功很多刚开始接触遥感数据的同事问起我都会把这段经验重新整理一遍。Himawari-8/9葵花8号、9号是日本气象厅的静止气象卫星定点在 140.7°E 赤道上空搭载 AHIAdvanced Himawari Imager成像仪一共 16 个波段全圆盘观测周期 10 分钟。对台风监测、强对流天气分析、火山灰扩散追踪、空气质量研究和海洋遥感来说它是亚太地区绕不开的数据源。这篇文章不打算翻译官方手册而是从你真的要把这些数据处理成能交付的东西这个角度出发把选数据源、读 HSD 格式、做定标、出真彩色图、批量处理这些环节里最关键的细节和坑讲清楚。如果你刚入门遥感或者已经被 HSD 文件折磨过几轮这篇文章应该能帮你省下大量试错时间。1. 先弄清葵花8/9的底细AHI的16个波段不是个个都该用于出图1.1 两颗卫星和它们的成像仪葵花8号2014年发射葵花9号2016年发射两者是同型号的在轨备份。目前实际业务运行的是9号但历史存档里大量数据来自8号。对做数据处理的人来说两者HSD文件结构、波段设置基本一致差别只在文件名的卫星编号H08/H09和个别标定系数的更新上。所以你可以把处理流程设计成对两者通用但文件命名和归档目录一定要分开不然后面做长序列分析时数据会乱成一锅粥。AHI 属于先进成像仪16 个波段覆盖可见光、近红外和热红外。和上一代 MTSAT 的5个通道相比信息量完全不在一个量级。也正因为通道多、频率高HSD 原始数据的体积和处理复杂度也跟着上来了——这是所有处理者要面对的第一个现实。1.2 波段光谱特性与用途速查处理这些数据的第一步不是写代码而是搞明白自己到底要哪个波段。下面是我自己整理的速查表做产品选波段时基本只看这一页波段中心波长(μm)星下点分辨率(km)主要用途适合出图时段B010.471蓝色可见光气溶胶、真彩色蓝通道白天B020.511绿色可见光真彩色绿通道白天B030.640.5红色可见光云识别、真彩色红通道白天B040.861植被与云近红外白天B051.62区分冰雪与云白天B062.32云粒子相态辅助白天B073.92火点、夜间云、海温全天B08-B106.2/6.9/7.32水汽通道全天B118.62云相态全天B129.62臭氧全天B1310.42红外窗区云顶温度全天B1411.22红外窗区海温全天B1512.42红外分裂窗全天B1613.32云顶高度全天注意 B03 是唯一一个 0.5 公里分辨率的波段单文件体积最大。如果只是做云图展示没必要每次都全量下载。而真彩色只依赖 B01/B02/B03夜间又完全用不上可见光要改用红外波段——这是不少新手踩的第一个坑晚上拿着一堆可见光数据最后只能对着黑屏发呆。1.3 先想清楚产品目标再谈处理在处理任何一批数据之前强烈建议先回答一个问题最终产出是什么是给业务汇报用的台风云图动画是用于反演海温的数据立方体还是做机器学习训练样本的真彩色贴图目标不同波段选择、定标级别、存储策略完全不同。举个例子如果你只要可见光真彩色图片那就只下载 B01/B02/B03 三个波段如果你要做云顶温度产品需要 B13甚至配合 B14/B15 做分裂窗如果是夜间监测优先 B07 和 B13 的组合。盲目把 16 个波段全部下载并解析磁盘会被快速塞满计算时间成倍增加最后大部分数据根本用不上。先定产品再定波段再定下载量这个顺序不要反。1.4 全圆盘10分钟一次意味着什么葵花8/9 每 10 分钟出一轮全圆盘观测。B03 单文件约 11000×11000×2 字节0.5km 分辨率B01/B02/B04 是 5500×5500×2其余 2km 波段是 2750×2750×2。一轮全波段全圆盘大约 600MB 上下一天 144 轮就是 85GB 量级。这个数字直接影响两件事一是下载必须做规划二是想长期保存某个波段的全历史数据存储成本要提前算。我的习惯是按需下载只保留高频使用的波段和区域原始 HSD 不长期囤积真正需要复算的时段再重新抓取。2. 数据从哪来下载渠道、HSD文件名和目录规划2.1 渠道对比与选用习惯葵花8/9 的数据分发渠道其实不少但每个渠道的定位和时效不一样。我用过的几个主要渠道如下渠道性质时效性适合场景JMA 官方分发正式数据需注册申请有延迟科研引用、业务结算、最终核验JAXA P-Tree 系统FTP/HTTP 公开归档准实时也有历史存档批量下载长时间序列高知大学气象信息网站HTTP 公开准实时快速拿当天数据做可视化NICT 实时分发网络HTTP 公开准实时快速查看、教学演示我个人的十字方针是批量回溯用 P-Tree快速出图用高知大学正式引用前从 JMA 官方渠道再核验一遍。数据源的解析逻辑都一样看图不挑渠道但发布、引用的正规材料必须挂可靠的官方来源这是行业基本要求。2.2 HSD文件名里的信息量每个 HSD 文件的名字已经把绝大部分元信息告诉你了。拿一个典型文件名拆解HS_H08_20190101_0000_B01_FLDK_R10_S0110.DAT片段含义HSHimawari Standard DataH08卫星编号H08 为葵花8号H09 为葵花9号20190101_0000观测时间UTC 时间 2019年1月1日00:00B01波段号FLDK观测区域Full Disk 全圆盘日本区域快速扫描文件常以 JP 开头如 JP01/JP08 等不同扫描模式编号R10空间分辨率R05 为 0.5kmR10 为 1kmR20 为 2kmS0110分段信息第01段共10段文件名里最容易出错的是时间。HSD 文件名用的是 UTC和北京时间差 8 小时和日本标准时间差 9 小时而且日本不实行夏令时换算关系全年固定。时区转换不能只做一次我的做法是所有脚本统一按 UTC 处理和归档只在最后出图标注时转成本地时间这样逻辑永远不会乱。2.3 下载与归档的目录习惯葵花数据量大如果没有目录规范三个月后你自己都找不着数据。我推荐这种布局himawari/ hsd/ 2019/ 01/ 01/ 0000/ HS_H08_20190101_0000_B01_FLDK_R10_S0110.DAT ...在根目录放一个manifest.csv记录已下载文件的名称、大小、本地路径和校验值。下载脚本每次启动先读 manifest已存在的文件跳过残缺的重新抓取。这样即使中途断网或磁盘满也能快速定位缺口。提示文件名时间是 UTC不是北京时间。所有自动化脚本请统一在 UTC 下工作最后出图时再转换能避免一大批看起来没毛病但就是不对的诡异问题。3. 手写HSD解析把二进制变成反射率和亮温3.1 HSD文件的基本结构HSD 每个波段单独一个文件整体结构可以拆成三层文件开头是 512 字节基本信息块记录卫星号、波段号、观测时刻、分段数、行列数等关键元数据之后是按块组织的数据区每个数据块前面有 128 字节的块头信息每个块内是按行排列的图像数据每一行由行号、16 位无符号整数像素、校验字节组成。理清这三层手写解析器并不难。难的是具体字节偏移和字段长度必须以 JMA 公布的 HSD 说明书为准不同版本可能存在差异。下面给一个示意性读取函数字段偏移写死之前一定要用你手上的真实文件验证一遍。3.2 Python读取HSD的骨架代码import numpy as np def read_hsd_full_disk(filename, ncol5500, nlin5500, seg_lines550, nseg10): data np.zeros((nlin, ncol), dtypenp.uint16) with open(filename, rb) as f: f.seek(512) # 跳过基本信息块 for seg in range(nseg): f.seek(128, 1) # 跳过本段块头 for row in range(seg_lines): f.seek(4, 1) # 跳过行号 raw f.read(ncol * 2) data[seg * seg_lines row, :] ( np.frombuffer(raw, dtypeu2) ) f.seek(4, 1) # 跳过行尾部校验 return data这个函数读完后data里就是原始 DN 值。验证解析是否正确有个土办法把数组做一次线性拉伸后直接存成 PNG如果能看到明显的云型结构说明行号偏移和每行字节数基本是对的如果图像全是斜条纹或噪声先回去查块头长度和行号字节数。3.3 可见光波段的定标HSD 里存的原始 DN 不是物理量。可见光和近红外波段需要通过增益和偏移转成反射率reflectance DN * gain offset这里的 gain 和 offset 在文件基本信息块里给出并且会随仪器状态修订。不要拿网上抄的固定值去套不同日期的文件否则时间序列上会出现系统性的跳变。做完反射率之后还要做太阳天顶角校正才能让不同时刻的观测可比reflectance_corrected reflectance / max(cos(solar_zenith), 0.05)太阳天顶角根据观测时刻、像素经纬度和太阳位置计算。不做这一步同一片云在早晨和中午的反射率差异巨大做序列分析时会出现假变化。分母加个 0.05 的下限是为了避免太阳在地平线附近时分母趋近于零导致数值爆炸。3.4 红外波段的亮温计算红外波段的 DN 同样先通过增益/偏移转成辐射值再用普朗克反函数转成亮温。核心公式是def radiance_to_bt(rad, wavenumber): c1 1.191042e-5 # W/(m2 sr cm-1)须与辐射值单位匹配 c2 1.4387752 # cm K return c2 * wavenumber / np.log(c1 * wavenumber**3 / rad 1)B13 中心波数约 947.95 cm⁻¹B14 约 891.34 cm⁻¹B15 约 801.61 cm⁻¹反演哪个通道就代入哪个中心波数。需要特别提醒辐射值的单位不同普朗克常数 c1 的数值也得跟着调整不然算出来的亮温可能会偏出几十K。这部分最容易出错我建议先找官方标记的样本数据验算一遍再进入批量。3.5 手写还是用现成库手写解析最大的价值在于理解链路出问题时知道该怀疑哪一环。但日常批量生产中我更推荐直接用现成框架也就是下一节要重点说的 Satpy。两者不冲突手写时积累的理解正好用来排查 Satpy 给出的奇怪结果。4. 提高效率的关键用Satpy直接出产品4.1 为什么选SatpySatpy 是 Pytroll 社区的开源遥感数据处理框架内置了 AHI HSD 格式的读取支持。它把读文件、定标、投影、合成这些环节都封装好了底层用 xarray 和 dask 做懒加载不会一次性把全圆盘数据灌进内存。对于没有时间手写解析器、又需要快速产出图像的场景Satpy 几乎是最优解。4.2 十五分钟出一张真彩色安装很直接pip install satpy出图代码比想象中短得多from glob import glob from satpy import Scene base /data/himawari/hsd/20190101/0000/ filenames glob(base /*B0[123]*FLDK*.DAT) scn Scene(filenamesfilenames, readerahi_hsd) scn.load([true_color]) scn.save_dataset(true_color, filenameh8_20190101_0000_tc.png)这段代码做的事包括识别 HSD 波段、读取基本信息、计算太阳位置、做反射率定标和太阳天顶角校正、执行真彩色合成、显示增强最后写 PNG。顺利的话十分钟到十五分钟就能出图。第一次跑通时你会觉得前面手写解析器的时间没白费因为遇到问题至少知道该去怀疑哪一环。4.3 真彩色不是三波段堆RGB真彩色合成并不是把 B03 当红、B02 当绿、B01 当蓝直接叠。原因有两个一是人眼对线性反射率的感知不是线性的直接叠加的结果暗部一团黑、亮部过曝二是大气分子散射尤其是蓝光会盖住地表细节让海面发灰、陆地发白。Satpy 的true_color合成器内置了大气订正和显示增强。想手动调优可以通过enhancement参数控制scn.save_dataset( true_color, filenameh8_tc_enhanced.png, enhancementgamma )伽马值一般取 2.0 左右但陆面、海面、晨昏边界的主观观感差异很大多试几个值才能找到顺眼的。这个环节没有绝对标准交付给不同业务方向口味也不一样。4.4 其他好用合成产品除了真彩色Satpy 里还有几个现成产品接业务时非常省事产品名说明典型用途natural_color自然色对冰雪更友好冬季云雪区分night_microphysics夜间微物理组合夜间云相态、雾区监测dust沙尘增强沙尘暴监测airmass气团分析产品锋面、急流分析ash火山灰增强火山灰监测这些产品不需要自己配波段公式scn.load([dust])就能用。应急监测场景里这套现成组合能帮你争取大量时间先出图再研究细节比从零开始搭公式靠谱得多。5. 区域裁剪、动画合成和地理叠加让数据真正能交付5.1 按经纬度裁剪全圆盘图很大但大多数时候只需要关注某片海域。Satpy 的crop方法可以直接按经纬度框选scn Scene(filenamesfilenames, readerahi_hsd) scn.load([true_color]) cropped scn.crop(ll_bbox(110, 10, 160, 45)) # 最小经度, 最小纬度, 最大经度, 最大纬度 cropped.save_dataset(true_color, typhoon_area.png)先裁剪再合成计算量会小很多。如果目标窗口特别小甚至可以只下载覆盖该范围的 HSD 分段从传输层就减少数据量。这个优化在带宽有限的环境下非常实用。5.2 生成动画的两条经验时间序列动画是台风复盘最常用的输出。我的流程是先把每个时次存成命名规范的 PNG再用 ffmpeg 合成视频ffmpeg -framerate 10 -pattern_type glob -i frames/*.png -c:v libx264 -pix_fmt yuv420p typhoon_2019.mp4这里有两个容易踩的坑。一是文件命名必须补零frame_1.png、frame_2.png……frame_10.png会被 glob 按字符串排成 1、10、2、3动画直接跳帧。正确命名是frame_0001.png这种。二是帧率决定视觉节奏10 分钟间隔的云图10fps 偏快5fps 更舒服具体看你想强调变化速度还是展示演变细节。5.3 叠加海岸线和地理标注单纯云图有时候不够直观叠上海岸线和经纬网格才方便定位。我通常分两步先把裁剪结果存成带坐标的 GeoTIFF再用 Cartopy 或 GIS 工具叠加岸线、城市名和台风路径。cropped.save_dataset(true_color, typhoon_area.tif, writergeotiff)拿到 GeoTIFF 后下面的脚本可以快速叠海岸线出图import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt import rasterio from rasterio.plot import show fig, ax plt.subplots( figsize(10, 10), subplot_kw{projection: ccrs.PlateCarree()} ) with rasterio.open(typhoon_area.tif) as src: show(src, axax, transformccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, lw0.6) ax.gridlines(draw_labelsTrue, linestyle--, alpha0.5) plt.savefig(typhoon_area_map.png, dpi150, bbox_inchestight)这一步不太复杂但不少做算法的同事会忽略导致交付的图看起来缺了地理感。加上海岸线之后整张产品的可读性完全不一样。6. 批量处理时的工程化问题和避坑清单6.1 磁盘、内存与下载策略前面算过全波段全天数据量约 85GB。如果只做真彩色动画只下载 B01/B02/B03 三个波段一天降到约 20GB能省四分之三的存储和下载时间。处理时优先用裁剪窗口减少内存占用0.5km 分辨率的 B03 整幅转成 float32 单波段接近 500MB多波段同时载入会让一般配置的机器直接卡死。利用 Satpy 的懒加载机制只在真正需要计算和保存时触发数据读取别上来就np.array()一把梭。6.2 时间、命名与维护窗口的坑一个很隐蔽的问题是时区。HSD 文件名的 UTC 时间在北京时间基础上减 8 小时和日本时间差 9 小时。我用的办法是所有文件名解析、归档目录、绘图标注统一用 UTC 字符串只在最后出图标题里转成本地时间。这样目录看起来晚了 8 小时但逻辑永远是自洽的不会出现某张图找不到对应帧的情况。另一个问题是卫星例行维护。葵花卫星会定期维护维护期间部分时次缺失。批量抓取前先检查当天文件数少于预期时先跳过不要拿残缺数据硬跑——否则动画里会突然出现一个黑屏帧汇报时解释起来很尴尬。6.3 常见症状与排查对照症状大概率原因处理方式真彩色全黑选了夜间时刻可见光无数据换红外产品或确认时间是否为 UTC图有规则横条纹分块/分行字节偏移算错校验行号字段对照官方说明修正 seek 长度颜色发红或发蓝波段顺序/合成公式不对检查使用的波段序列是否正确B03/B02/B01动画时间乱跳文件名按字符串排序、没补零用补零命名或按解析出的 UTC 时间排序下载后文件数量不足网络抖动或上游维护做文件数自检缺失时重新抓取6.4 数据引用与合规使用葵花数据时交付报告和出版物里要按 JMA 要求标注数据来源标准写法是注明 Original data from Japan Meteorological Agency。这是行业里的基本礼貌也是科研伦理要求。归档时保留 HSD 文件的原始文件名不要擅自改名方便回溯和复现。最后再分享一个我从踩坑中换来的习惯任何环节的脚本都先拿单个文件、单个时次跑通并保存输出确认结果完全符合预期后再放进批处理循环。卫星数据处理的链路很长如果一次性跑几千个文件出了问题你根本分不清是定标问题、投影问题还是文件名解析问题。先小后大、先单后批看起来慢实际操作下来反而是最快到终点的路。这套流程我现在无论处理葵花还是其他静止卫星数据都一直在用。