简介这是一套使用Python读取新一代天气雷达基数据的程序包面向气象科研人员、雷达数据处理工程师及有Python基础的学习者覆盖雷达基数据导入、清洗、质量控制、格式转换、分析与可视化等全流程。资源共6个文件包含5个Python脚本和1个说明文档压缩包仅8KB代码规模精简、便于阅读理解。脚本按功能划分清晰覆盖数据导入、质量控制、标准化转换、识别分类等处理环节附带说明文档提供使用指南与功能简介辅助用户快速上手。已有103人学习下载该资源。读者可由此掌握新一代天气雷达基数据的标准读取思路理解数据质量控制、标准化格式转换、雷达数据识别分类等关键环节的实现方式为后续开展雷达资料分析、算法研究或二次开发提供一套可直接参考复用的Python代码框架。1. 新一代天气雷达基数据为什么非得自己写Python程序来读某天某省台的同行传过来一个雷达基数据文件后缀是 .bin说帮我画一张反射率图。我接过文件用 numpy 直接读出来的全是乱码偶尔能看见 CINRAD 几个 ASCII 字符像是有规律又抓不住。读取新一代天气雷达基数据Python程序做的正是这件事把 SA/CB 型新一代天气雷达下发的基数据二进制流解析成 numpy 数组、物理量场供后续画图、拼图、定量估测降水使用。为什么必须自己写因为基数据不是通用格式市面上的开源雷达库大多面向国外格式落地到国内站点的基数据上基本是黑匣子。这篇文章写给要做回波定量分析的开发、要做短临算法的学生、以及想彻底搞懂雷达基数据底层结构的人读完能照着复现一套完整读取流程。2. 基数据的文件结构文件头、仰角块与径向记录的读取顺序2.1 基数据和产品数据是两回事为什么必须读基数据新一代天气雷达的输出分两级。一级是基数据也就是雷达体扫过程中逐脉冲累积、平均、退模糊后得到的原始观测包含反射率、径向速度、速度谱宽三个核心量另一级是产品是基数据经过插值、坐标转换、组合反射率等算法加工后的图像或网格文件。产品看起来直观但坏了关键的一点算法是雷达厂商封死的你想改阈值、改插值半径、做自定义拼图产品数据给不了你自由度。做短临预报和定量降水的人最终都要回到基数据上自己动手。反射率可以算雨强速度场可以反演风场谱宽可以判断湍流这三个量在基数据里全是原始采样值时间分辨率能到 6 分钟一个体扫空间上每条径向有好几百个距离门。要拿到这些只能解析基数据文件。常见做法是用 Python 的 struct 模块按字节解包配合 numpy 组织成多维数组这也是本文章标题所指程序的核心工作。2.2 文件布局一条记录一条径向靠识别码切分我拆过好几个站点导出的基数据文件整体布局基本一致。文件开头是文件头区记录雷达站号、扫描时间、体扫模式、仰角层数、每个仰角的径向数接着是若干仰角数据块每个块内按顺序存放该仰角的所有径向记录每一条径向记录又由记录头和样本数据组成。记录头里写的是方位角、仰角、径向样本数、数据起始位置数据段就是反射率、速度、谱宽的编码值。雷达型号不同字段偏移量会差几个字节但文件的 ASCII 标识是稳定的。几乎所有国产新一代天气雷达基数据文件头里都有一段 ASCII 文本形如 CINRAD后面的字节串里能直接读出站号和雷达型号。我一般先读整个文件到内存再用 find 定位这个标识而不是硬编码偏移量。下面是常见文件布局的对照表区域内容读取要点文件头站号、型号、时间、体扫模式、仰角层数用 ASCII 标识定位仰角索引块每个仰角的径向数、起始字节偏移动态读取不要写死径向记录头方位角、仰角、样本数、数据长度小端字节序结构体对齐径向数据段反射率 1 字节、速度 2 字节、谱宽 2 字节按样本数切片即可2.3 为什么用 Python 而不是 C 或 Java解析二进制这件事C 做起来性能最优但气象算法迭代快用 C 改一版解码逻辑要重新编译调试成本高。Java 的字节处理不差但生态里缺少 numpy 这种多维数组基础设施后面做径向数据重组、极坐标转直角坐标都得自己造轮子。Python 的 struct 加 numpy 是业界最常见的组合struct 负责按格式字符串解包固定长度的二进制块numpy 负责把解出来的径向数组一次性堆成二维网格。这里有一个重要习惯解析基数据不要逐字节 read而是一次性把整个文件读进内存再用切片取数据段。基数据单个文件通常几十 MB一次读完完全没压力还能避免频繁 IO 带来的文件句柄问题。字节序方面基数据默认小端struct.unpack 的格式串要记得加前缀比如H表示无符号短整型。真遇到个别大端导出的文件解出来方位角会大得离谱这种情况在避坑章节细说。3. 用Python解析径向数据从字节流到反射率、速度、谱宽场3.1 先定位文件头别上来就硬算偏移量解析基数据的第一个动作是把文件读成 bytes然后找 ASCII 标识。我踩过最深的坑就是拿某一站的基数据反推了偏移量代码在本站跑得完美换一个站直接乱码。后来改成扫描式定位问题再没出现过。先看定位函数import numpy as np def read_into_bytes(filepath): with open(filepath, rb) as f: data f.read() return data def find_cinrad_header(data: bytes) - int: # 扫描文件中的 CINRAD 标识返回标识的起始字节偏移 idx data.find(bCINRAD) if idx 0: raise ValueError(未找到 CINRAD 标识请确认是新一代天气雷达基数据文件) return idx逻辑说明read_into_bytes一次性载入全部字节避免多次 IOfind_cinrad_header用 bytes.find 扫描标识。文件头不一定在文件最开头因为部分雷达站会先写一段固定字长的使用说明或填充区所以扫描比固定偏移稳。参数说明data.find(bCINRAD)返回第一个匹配位置找不到返回 -1。如果站点导出的基数据标识是 CINRAD/SA 之类find 依然能命中前缀。找到标识后从标识后连续取 6 字节可读出站号再取 4 字节读出雷达型号编码格式多为 ASCII注意去掉末尾的\x00填充。3.2 解析一条径向记录头与数据段的拆包定位到文件头之后真正的工作在径向层面。每一条径向是一个独立的数据单元我从实际解析经验出发建议先写一个只解析单条径向的函数验证通过后再套循环。下面是我常用的解析函数import struct def parse_radial(data: bytes, offset: int): # 记录头前 4 字节为记录长度之后是方位角、仰角、距离门数 header data[offset:offset 32] radial_len, az_index, el_index, gate_count struct.unpack(HHHH, header[0:8]) az struct.unpack(H, header[8:10])[0] / 4096.0 * 360.0 # 方位角 el struct.unpack(H, header[10:12])[0] / 4096.0 * 90.0 # 仰角 # 反射率段每距离门 1 字节 dbz_start offset 32 dbz_raw np.frombuffer(data, dtypenp.uint8, countgate_count, offsetdbz_start) # 速度段每距离门 2 字节紧接反射率之后 vel_start dbz_start gate_count vel_raw np.frombuffer(data, dtypenp.uint16, countgate_count, offsetvel_start) # 谱宽段每距离门 2 字节紧接速度之后 sw_start vel_start gate_count * 2 sw_raw np.frombuffer(data, dtypenp.uint16, countgate_count, offsetsw_start) return { azimuth: az, elevation: el, dbz: dbz_raw, velocity: vel_raw, spectrum_width: sw_raw, }逻辑说明struct.unpack(HHHH, header[0:8])一次解出四个无符号短整型分别是径向记录总长、方位角索引、仰角索引、距离门数。方位角索引和仰角索引不是角度本身而是 0 到 4096 的整数映射换算时除以 4096 再乘对应量程这是基数据常见的定点表示法。数据段用np.frombuffer直接映射内存不产生 Python 列表速度上有数量级优势。参数说明dtypenp.uint8对应反射率一个距离门一个字节dtypenp.uint16对应速度和谱宽各占两字节。offset参数告诉 frombuffer 从哪个字节开始读取。注意vel_start的计算依赖gate_count所以记录头解析出错时速度段会整体错位后期画图出现风场断裂排查时先检查 gate_count。3.3 组装单仰角径向数据把多条径向拼成二维场解析出单条径向还不够要得到一个仰角层的完整扫描需要把该仰角下所有径向按方位角拼成二维数组。这里要理解雷达扫描的极坐标特性每一条径向是一根射线距离门从雷达中心向外延伸不同射线方位角不同。组装函数如下def assemble_sweep(radials): # radials 是 parse_radial 返回的字典列表 n_radials len(radials) gate_count radials[0][dbz].shape[0] dbz_field np.full((n_radials, gate_count), np.nan, dtypenp.float32) az_array np.zeros(n_radials, dtypenp.float32) for i, r in enumerate(radials): dbz_field[i, :] r[dbz] az_array[i] r[azimuth] return az_array, dbz_field逻辑说明先用 np.full 初始化一个全 NaN 的二维数组避免数据不足时出现 0 值的假象。逐径向把字节数组复制进对应的行方位角单独存一列数组后续做极坐标绘图时用作坐标轴。这里没有按方位角排序因为基数据本身是按扫描顺序存储的先按方位角大小排序会让绘图更规整但会影响与原始径向的对应关系我一般在绘图时排序在数据层保持原始顺序。参数说明np.nan是缺测值的统一标记后续做反射率等值线、统计面积时会自动跳过dtypenp.float32用于精度和内存的折中一个 1000x400 的场约 4MB。az_array的角度单位是度范围 0 到 360绘图前如需从正北开始要做偏移换算这也是新手最容易忽略的坐标细节。4. 一次读完整个体扫多仰角组织、批量文件与物理量解码4.1 体扫数据的组织方式从单仰角到三维结构单仰角组装完成之后可以进一步组织体扫数据。体扫包含多个仰角层每层径向数相同但距离门数可能不同。组织体扫时我用一个字典结构键是仰角度数值是组装好的二维场和方位角数组这样上层代码可以直接按仰角取值def parse_volume(filepath): data read_into_bytes(filepath) idx find_cinrad_header(data) # 从文件头读取仰角层数通常位于标识后第 18 字节按小端读 n_sweeps struct.unpack(H, data[idx 18:idx 20])[0] elev_deg struct.unpack(HH, data[idx 20:idx 24]) # 示例具体以实际站点为准 sweeps {} offset idx 256 # 头部预留区长度不同版本有差异 for _ in range(n_sweeps): radials [] for _ in range(radials_per_sweep): r parse_radial(data, offset) offset r[dbz].shape[0] * 1 r[velocity].shape[0] * 2 r[spectrum_width].shape[0] * 2 32 radials.append(r) az, dbz_field assemble_sweep(radials) # 以该仰角数据中的第一个仰角值为键 key round(radials[0][elevation], 1) sweeps[key] {azimuth: az, dbz: dbz_field} return sweeps逻辑说明parse_volume先定位标识再从文件头读层数之后进入双循环。外层循环走仰角内层循环走径向每次循环后手动推进 offset推进量等于记录头长度加三个数据段长度之和。这里的n_sweeps和radials_per_sweep都应从文件头动态读取实际项目中我会把这两个值解析出来而不是用常量。参数说明文件头预留区长度 256 只是示例值不同雷达站导出的文件在头部预留区长度上有较大差异。如果发现解出来的仰角值或者层数异常优先检查这个预留区长度。key round(radials[0][elevation], 1)用仰角值做字典键可以避免浮点误差造成的键不匹配。我还会额外保存原始字节数据到字典里方便排查时逐字节比对。4.2 把索引值变成物理量反射率、速度、谱宽的解码基数据文件里存的不是直接的物理量而是索引值。反射率、速度、谱宽各自有一套映射关系读取程序里最关键的一步是把这些索引还原成气象上可用的物理量。反射率通常直接查表速度和谱宽则一般通过索引与缩放系数相乘得到实际值。下面是我常用的解码函数def decode_dbz(dbz_raw, offset_db0.0): # 反射率索引值映射到 dBZ常见映射为 索引/2 - 32 # 索引 0 和 255 通常是无效值映射为 NaN valid (dbz_raw 0) (dbz_raw 255) dbz np.full_like(dbz_raw, np.nan, dtypenp.float32) dbz[valid] dbz_raw[valid] / 2.0 - 32.0 offset_db return dbz def decode_velocity(vel_raw, vel_res0.0625): # 速度索引值乘以分辨率得到 m/s索引 0 或 65535 表示缺测 valid (vel_raw 0) (vel_raw 65535) vel np.full_like(vel_raw, np.nan, dtypenp.float32) vel[valid] (vel_raw[valid] - 32768.0) * vel_res return vel逻辑说明decode_dbz里把 0 到 255 的索引值除以 2 再减 32得到一个常见的 dBZ 映射区间索引 0 和 255 视为无效值丢弃。这里的offset_db参数可以用来做雷达标定偏差修正。decode_velocity的写法是先把索引减去 32768 变成带符号数再乘以速度分辨率这样既支持负速度也支持正速度符合径向速度的定义。参数说明vel_res0.0625是常用速度分辨率但 S 波段和 C 波段雷达有时候不一样。这里有玄学成分同一个站不同时期升级后分辨率都可能变拿到数据后应该先用已知的台风或强对流过程对比一下最大速度值是否合理。np.full_like保留了输入数组形状并且对无效索引自动置 NaN这一步不做的话后续画图会出现大面积 0 值误导。4.3 批量文件与时间序列一次解析一个时次的全部文件实际业务中往往收到一个时次的多个体扫文件或者一个连续过程的全部文件。批量读取时要按文件名排序保持时间顺序。下面是批量处理函数import glob import os def batch_parse(base_dir, file_pattern*.bin): file_list sorted(glob.glob(os.path.join(base_dir, file_pattern))) volume_list [] for fp in file_list: vol parse_volume(fp) # 从文件头或文件名中提取时间戳存到结果里 vol[time] os.path.basename(fp)[:14] # 常见命名YYYYMMDDHHMMSS... volume_list.append(vol) return volume_list逻辑说明batch_parse的关键有三点。第一glob模式匹配后必须 sorted因为文件系统返回顺序不保证按时间排列第二时间戳从文件名前 14 位截取前提是站点导出的文件按 YYYYMMDDHHMMSS 命名这个规则在各省略有差别读取前先人工看一眼文件列表第三每解析完一个体扫把结果存进列表不要全部累积到内存再处理文件多时内存容易爆。参数说明file_pattern*.bin可以按实际后缀调整有的站点导出是.dat或没有后缀。vol[time]我建议存成字符串而不是 datetime 对象避免解析失败拖慢整个批量过程。处理完所有文件后可以把volume_list用 pickle 或 np.save 存下来下次调试算法时直接加载省去重复解析的时间。5. 解析基数据的五个典型坑现象、原因、对策5.1 换一个站点的基数据解出来的方位角乱跳现象程序在 A 站数据上跑得正常换到 B 站的同型号文件后方位角出现大于 360 度或者负数的异常值反射率场形如乱码。原因不同站点导出的基数据文件头预留区长度不同我最初拿 A 站文件反推的固定偏移量在 B 站直接错位导致径向记录头里的方位角字段读到了数据段的字节。这不是解析逻辑的问题而是头部定位方式的缺陷。解决改成 3.1 节里的扫描式定位用data.find(bCINRAD)定位标识后再逐字段解析。如果标识后字段仍然错位就在标识后每隔 4 字节尝试解一次站号和仰角层数校验通过再继续。这套方法看起来笨但对付各站点格式差异是最稳的。5.2 速度场全是毛刺风羽乱飞现象速度场画出来后回波区内部出现密集的正负交替斑块看起来像噪声与当天实况风场明显不符。原因速度索引到 m/s 的换算关系不对。不同雷达型号、不同体扫模式下的最大可测速度不同速度分辨率也不固定直接套用默认分辨率会导致速度值被放大或缩小若干倍。我是拿一次已知的层状云降水过程做过对比才发现手里的分辨率参数过时了。解决先用本站的历史基数据比对反射率与速度的对应关系找到合理的最大速度值反推分辨率更稳的做法是找雷达站标定文档里的速度查找表直接把索引映射成数组查表。代码里要把换算参数提出来做成配置项方便一个站点一套参数。5.3 仰角层数和文件头写的不一致现象程序读完第一层后循环进入第二层时索引越界或者第二层仰角值变成 0。原因基数据的体扫模式决定仰角层数但有些文件头里的层数是预置值实际数据块里只写了几层就收尾了。比如文件头写了 9 层但数据仅 5 层有效后 4 层是空记录。解决解析时先用循环找记录头里的仰角字段如果仰角值小于 0 或大于 90判定为无效径向停止该仰角片段的读取再校验层数与文件头声明的差异以实际读取到的有效层为准。追加一层保护逻辑宁可少一层不可越界崩溃。5.4 反射率色标整体偏亮或偏暗现象反射率图整体比产品图偏亮原本 30 dBZ 的回波区域在图上是 40 dBZ整场偏移。原因反射率索引到 dBZ 的映射表不同。部分雷达使用非线性映射索引值并非简单的线性对应直接用索引除以 2 减 32 会把中间段的值整体抬高。这个坑在弱回波区表现不明显强回波区偏差很大。解决让用户提供一份站点反射率查找表把每个索引值对应的 dBZ 写进配置文件程序启动时加载并查表没有查找表时用线性映射但保留offset_db参数人工修正。还有一个土办法把解析结果与本站产品图上的回波顶高做对比反推偏差量。5.5 大文件解析慢一个体扫要十几秒现象单个体扫文件解析完耗时明显偏高批量处理一个过程十几个文件跑到最后越来越慢怀疑内存泄漏。原因逐径向使用struct.unpack会创建大量临时对象组装径向时反复调用np.array也有开销。另一个常见原因是文件被反复 seek 读取而不是一次性读入内存。解决一次读入全部字节后用np.frombuffer加 offset 直接切片避免逐字节解包循环内避免创建 Python 列表装字节数组直接用 numpy 索引切片。实测这样能把单文件解析压到 1 秒内。如果数据量再大用 memoryview 包一层切片减少内存复制。6. 读取之后怎么验证和提速把基数据转成常用格式再干活6.1 画一张 PPI 图验证解析正确性解析程序写完先别急着做算法先画一张 PPI平面位置显示图验证场的空间形态是否正确。用 matplotlib 的 pcolormesh 配合极坐标投影是最直接的验证手段import matplotlib.pyplot as plt def plot_ppi(az, dbz_field, max_range_km230): # az 为方位角数组dbz_field 为反射率二维场 r, theta np.meshgrid(np.linspace(0, max_range_km, dbz_field.shape[1]), np.deg2rad(az)) fig, ax plt.subplots(subplot_kw{projection: polar}) mesh ax.pcolormesh(theta, r, dbz_field, cmapjet, vmin0, vmax60) ax.set_theta_zero_location(N) ax.set_theta_direction(-1) plt.colorbar(mesh) plt.show()这段代码里np.deg2rad(az)把方位角转成弧度set_theta_zero_location(N)把 0 度放到正北set_theta_direction(-1)让角度顺时针增加符合雷达扫描习惯。如果画出来的回波是半个圆弧或者位置偏移先查方位角换算再查距离门的长度单位。6.2 序列化为 npy 和 xarray避免重复解析基数据解析的时间成本虽然能压到 1 秒内但批量调试算法时反复解析仍然痛苦。我一般解析完直接存成压缩 npy 或 xarray后续代码直接加载。单文件体扫存 npy 的方式np.savez_compressed(vol_20230712.npz, **sweeps)加载时用np.load后按字典取键即可。如果后续要做坐标变换或与数值模式数据对照推荐用 xarray 组织把方位角、距离门、仰角都标成维度坐标插值和切片会顺手很多。存储后的加载速度比重新解析快一个量级调试迭代的体验好很多。6.3 验证基数据质量后再谈业务应用读取程序稳定后下一步拿它做业务之前先跑一遍数据质量体检统计每个仰角层的有效径向数、缺测比例、速度退模糊后的最大速度是否合理再与雷达站发布的对应产品图对比确认方位角和距离门没有系统偏差。这步做完后续的拼图、估测降水、风场反演才有可信度。这些年我见过太多解析程序跑通就跑业务最后在极端个例上翻车的案例。解析基数据这件事本质上是在和二进制格式做斗争格式版本一变之前的偏移量全部作废。我现在的习惯是解析层和业务层严格分离解析层只输出原始索引和标准物理量业务层永远只吃解析后的结果每一份新站点的基数据先跑一遍校验脚本再入库。这个习惯帮我少踩了很多坑希望帮到你。本文还有配套的精品资源点击获取