Python解析CINRAD雷达基数据：从二进制解码到PPI/RHI绘图实战

简介：基于Python的CINRAD雷达数据读取与绘图源码，是一套面向气象数据分析师、科研人员及高校学生的完整雷达数据处理方案。它解决CINRAD基数据读取难、可视化流程繁琐的问题，支持批量读取与交互式操作，可绘制PPI、RHI及多种产品图，并配套REF、VEL、ZDR等色标文件，适配气象业务分析与教学场景。资源共82个文件，以28个Python脚本为核心，涵盖数据解析、坐标投影、可视化与常用计算模块；另有14个颜色映射文件、6个XML配置、3个spec格式说明，以及含城市/省界的地理信息与雷达数据文件，压缩包约38MB。项目还提供示例脚本、Jupyter Notebook演示和中文说明文档，便于快速上手。目前已有936人学习此资源，适合需要掌握雷达数据读取与可视化、希望复用代码搭建分析流程的气象从业者与学习者。 干气象数据处理这行，应该都见过CINRAD雷达的基数据文件，这其实就是比较常见的雷达数据承载形式。中国新一代天气雷达网（CINRAD）在业务运行中持续产出这种二进制格式的观测数据，里面存的是反射率、径向速度、速度谱宽这些基础物理量，但文件本身直接用普通编辑器打开，就是满屏乱码。做科研和业务开发时，我们需要把这些数据读取出来，再绘制成常用的PPI（平面位置显示）图、RHI（距离高度显示）图，甚至叠加地理底图做分析。我就基于Python从零写了一套读取与绘图的源码，不依赖专门的闭源工具，也不调用现成的第三方CINRAD读取库，纯标准库加科学计算组件把二进制解码做完整。

这套源码适合的人很具体：气象专业学生做毕业设计、搞雷达资料同化或短临预报的科研人员，以及想理解雷达基数据存储结构的工程师。它解决两个实际痛点，一是把雷达二进制格式读干净，二是出图可控，从色标到显示范围都能自己定。下面我把整个设计思路和实现过程拆开讲，包括格式解析、物理量换算、画图逻辑和踩坑记录，希望能给正在折腾雷达数据的同学省点时间。

1. 需求拆解与技术选型

1.1 CINRAD基数据文件到底存了什么

CINRAD/SA雷达（多普勒天气雷达）的基数据文件是典型的二进制记录结构，一套完整时次的数据大致包括文件级头记录、径向数据块和尾记录三个部分。文件头记录里包含站点参数和雷达参数，比如站点编号、雷达经纬度、天线类型、发射机频率、脉冲宽度等；径向数据块是文件主体，每条径向存储仰角、方位角、距离库数，以及若干物理量的原始数值；尾记录负责数据完整性标记。里面最关键的是“每根径向的仰角、方位角、距离库序列”，有了这些，才能知道一根雷达射线从哪个角度发射、回波在哪些距离上出现，最终还原整个扫描过程。

为什么选Python来做这个事，背后是有实际考量的。雷达数据解码本质上是按固定字节偏移量去取值，Python的struct模块做二进制解析非常顺手；numpy又能把一根根径向数据堆叠成二维数组，后续画图直接用matplotlib处理。对比用C++或Fortran，Python开发周期短很多，对科研场景尤其友好，而且兼容性广泛，从Windows笔记本到Linux服务器都能跑。如果以后想把算法移植到业务平台，也可以用cython或nuitka打包，并不会有太高的迁移成本。

1.2 模块划分与整体设计思路

我设计源码的时候没有把所有逻辑堆在一个文件里，而是拆成三层。第一层是解码层，负责打开二进制文件，逐条读取记录，解析出物理量数组；第二层是数据层，把解码结果封装成体扫对象，包含仰角、方位角、距离数组、反射率数组等；第三层是可视化层，把体扫对象转换成图形，支持PPI、RHI等模式，负责插值、色标和地理叠加。

这样分层的理由很简单：解码逻辑和绘图逻辑应当彻底解耦。解码接口保持稳定，后面换绘图引擎或调整出图样式，完全不用动解码层。扩展也很方便，比如增加读取雨量计数据叠加、导出NetCDF，都可以在各层内独立完成。对这个项目来说，还有个容易被忽略但很关键的选择——不直接用现成的第三方CINRAD读取库。第三方库虽然能快速出结果，但黑盒特性太强，很难理解字节偏移和物理量换算的细节。自己实现解码，反而能把格式吃透，日后遇到非标准文件也有能力处理，这在科研中其实比“能用”重要得多。

2. CINRAD核心格式解析与解码实现

2.1 二进制结构与字节序

CINRAD基数据文件按记录组织，每一条记录往往以长度标记开头、结尾，这种结构类似很多通信协议里的TLV（类型-长度-值）设计。文件头部记录里包含站点编号、雷达经纬度、天线参数等信息，这些字段大多是小端字节序存储。解析时，Python最常用的组合是open()函数配合struct.unpack()按格式串解包。比如读取站点编号，可以先读取记录总长度（4字节小端），再按偏移量取站点字段。

import struct with open('RADA_CHN_YBZ_CINRAD_SA_CREF_20230720_000000.bin', 'rb') as f: # 先读取记录总长度（4字节，小端） rec_len = struct.unpack('<I', f.read(4))[0] # 站点编号偏移量需要参考格式文档 f.seek(6) station_id = struct.unpack('<4s', f.read(4))[0].decode('ascii', errors='ignore')

这里需要注意，格式化串里的小于号（<）表示小端，千万别漏掉。之前有同事拿C++程序移植到Python时直接把'<I'写成'I'，结果读出来的记录长度全部高出几万倍，定位问题就花了一个下午。另外，不同型号雷达（SA、CB、CC）文件头的字段偏移会有差异，务必以对应格式说明为准。解码的核心原则是：先拿到记录的总长度，再按字段偏移量逐一取读，最后用总长度校验是否读完整。很多解析乱码问题都是因为偏移量算错一位导致后续字段全部错位，整个文件数据就全乱了。

2.2 径向数据解包与物理量还原

径向数据块是体扫的主体，一个体扫文件里包含大量径向，每条径向都由仰角、方位角、距离库数和数据数组组成。解码之后要做的一件重要事情叫物理量换算，因为雷达原始存储的数值并不是直接可用的物理量，而是经过缩放的数据。常见换算关系如下：反射率因子实际值等于原始值除以10再减66，单位是dBZ；径向速度实际值等于原始值除以10再减127.5，单位是m/s；速度谱宽实际值等于原始值除以10再减63.5，单位是m/s。

换算时要同时处理无效值。雷达在无回波区域会存成特定无效标记，比如0或65535，如果不加处理直接带入换算，画出来的图会到处是噪点，而且会污染插值结果。我写了个读取径向数据的函数，伪代码结构如下：

import numpy as np def read_radial(f, gate_count): azimuth = struct.unpack('<H', f.read(2))[0] / 100.0 elevation = struct.unpack('<H', f.read(2))[0] / 100.0 raw_ref = struct.unpack(f'<{gate_count}H', f.read(gate_count * 2)) ref = np.array(raw_ref, dtype=np.float32) # 把无效值替换为nan，避免后续画图时被当作有效回波 ref = np.where((ref > 0) & (ref < 40000), ref, np.nan) ref = ref / 10.0 - 66.0 return azimuth, elevation, ref

实际格式里字段顺序可能不同，需要对照格式文档调整。有一个容易踩坑的细节：不同雷达型号对无效值定义不同，有的存0、有的存65535，所以我会在解码层做一个可配置的无效值范围，而不是写死单一阈值。这个设计后期帮我省了不少事，因为换了个省份的雷达文件试跑时，基本不用改代码。

2.3 从径向数据到三维体扫

把逐条径向数据组织成体扫对象时，我按仰角分组。每个仰角层包含若干方位角上的径向，距离库数量一般固定。这样自然得到三维数组：仰角乘以方位角乘以距离库。一个体扫的典型参数是9个体扫仰角、360根方位角每1度一根、每根径向460个距离库，存储成numpy三维数组后，后续绘图和统计都非常高效。

这部分还需要计算每个距离库对应的物理距离。距离库的长度主要由脉冲宽度决定，常见短脉冲的距离库长250米，反射率扫描一般采用1公里量级的库长。计算方式很简单，距离等于距离库序号乘以库长，但在做经纬度投影时还要结合雷达站的经纬度和每个仰角、方位角的波束传播路径。业务上多数时候不需要精细的波束弯曲修正，但做定量估测降水时，这个修正还是有必要做的，涉及到标准大气折射条件下的4/3地球半径近似。

3. 绘图模块设计与可视化实现

3.1 PPI图的两种实现思路

PPI图是把某个仰角层上的雷达观测，按方位角和距离展开到水平平面。初学时最容易想到的做法是，把每个距离库当作直角坐标上的点，直接用scatter散点画。但这样画出来的图密度不均匀，远处稀疏、近处密集，视觉效果很差，也不方便和地理信息叠加。我采用的方法是网格化插值：先把极坐标下方位角和距离的每个点转换到直角坐标，然后构建规则网格，再用最近邻或线性插值把散点数据映射到网格上，最后用pcolormesh或imshow填充颜色。

from scipy.interpolate import griddata import matplotlib.pyplot as plt # thetas: (num_azimuth, 1) 方位角，r: (num_gates, ) 距离库 lon_grid = r * np.cos(np.deg2rad(thetas))[:, np.newaxis] lat_grid = r * np.sin(np.deg2rad(thetas))[:, np.newaxis] points = np.column_stack([lon_grid.ravel(), lat_grid.ravel()]) values = ref_data.ravel() xi, yi = np.meshgrid(np.linspace(-r.max(), r.max(), 500), np.linspace(-r.max(), r.max(), 500)) zi = griddata(points, values, (xi, yi), method='linear') zi = np.ma.masked_invalid(zi) fig, ax = plt.subplots(figsize=(8, 6)) pm = ax.pcolormesh(xi, yi, zi, cmap='reflectivity', vmin=0, vmax=70) plt.colorbar(pm, ax=ax, label='dBZ') ax.set_aspect('equal')

这段代码在我的机器上能直接跑通，但有几个细节容易坑人。角度转弧度必须单位统一，我见过有人忘记把度数转弧度，画出来的图就像被人切掉了一块；网格范围务必以最大探测距离为边界，否则空白处会被插值填出很多假回波；插值方法选linear时，数据稀疏区容易出现边缘锯齿，这时候改成nearest会好一点，虽然颜色块感重一些。

3.2 色标设计与业务对齐

气象雷达图的色标不是随便拿个jet就完事了。业务上常用的反射率色标是蓝绿黄橙红紫渐次增强，速度图则是红绿两色区分正负径向速度，谱宽图偏蓝色系。为了不每次画图都重新调色，我把色标单独封装成一个模块，定义了反射率色标从-10 dBZ到70 dBZ、速度色标从-27 m/s到27 m/s、谱宽色标从0到10 m/s三套标准色标。

用matplotlib实现时，最直接的方法是自定义LinearSegmentedColormap。以反射率为例：

from matplotlib.colors import LinearSegmentedColormap colors = ['#00ffff', '#00ccff', '#0099ff', '#0066ff', '#0033cc', '#00ff00', '#00cc00', '#009900', '#ffff00', '#ffcc00', '#ff9900', '#ff6600', '#ff3300', '#ff0000', '#cc0000'] ref_cmap = LinearSegmentedColormap.from_list('reflectivity', colors)

这样做的收益很实际：和气象业务平台发布的图基本一致，审稿人和业务人员看起来不别扭，也省掉每次出图都反复解释颜色的麻烦。速度图尤其要注意，过零色一定要用浅灰色或者白色，不然正负速度分界看不清。另外，matplotlib 3.x版本中，imshow的vmin/vmax改成vmin/vmax后依然生效，但最好显式设置clim，这样在批量出图时能保持一致的颜色映射范围，防止横向对比时出现“同样的回波强度却用不同颜色”的尴尬情况。

3.3 RHI剖面与多仰角拼接

除了PPI，RHI在分析回波垂直结构时非常常用。RHI图沿一条径向切面，横轴是距离，纵轴是高度。高度计算需要考虑标准大气折射率条件下的波束高度公式，采用有效地球半径4/3近似，公式为h等于sqrt(r平方加kRe括号平方加2r kRe sin(elev))再减去kRe，其中r是斜距，elev是仰角，Re是地球半径，k取4/3。这个公式并不复杂，但很多初学者容易把斜距当成水平距离直接画，导致图上的回波被压缩在低层，看起来全是地形遮挡效果，实际上只是坐标算错了。

多仰角拼接则是把多个仰角层数据显示在一张图里，用于合成反射率。实际业务中最常用的Composite Reflectivity，就是把所有仰角层每个距离库取最大值，达到最大覆盖范围的展示效果。实现时，用numpy沿仰角维度做nanmax压缩，一行代码就能搞定。要注意的是，如果某个距离库在所有仰角层都是无效值，nanmax会返回nan，画图时需要用mask过滤，否则会出现一整片纯色块，误导读图者以为那里有强回波。

4. 源码架构与关键实现细节

4.1 目录结构与运行方式

我给出一套可直接落地的目录结构，方便直接对照参考。这里没有把类写得很复杂，全部用函数加轻量封装，适合阅读和二次开发。目录大致分成五个部分：解码层cinrad_parser.py、数据层radar_data.py、可视化层plot_radar.py、命令行入口main.py，以及data和output两个文件夹分别存原始文件与输出图。

cinrad_tools/ ├── cinrad_parser.py # 解码层：二进制读取与物理量换算 ├── radar_data.py # 数据层：体扫对象、径向分组 ├── plot_radar.py # 可视化层：PPI、RHI、色标 ├── main.py # 命令行入口 ├── data/ # 存放雷达基数据文件 └── output/ # 输出图像

入口main.py用argparse接收文件路径、产品类型、仰角序号等参数。运行方式很简单，在命令行敲：

python main.py -f data/RADA_CHN_YBZ_CINRAD_SA_CREF_20230720_000000.bin -p ppi -e 2 -o output/

这样设计的好处是便于批处理，写一个for循环就能对一整个时次的文件夹批量出图，不用每次手动改代码。我在做一次强对流过程分析时，就用这个入口一口气处理了48个时次，每个时次画三张图，耗时十几分钟，基本满足科研需求。

4.2 核心函数设计

解码层的核心是parse_volume函数，按记录号顺序依次读取和校验。我会把记录总长度放在循环开头和结尾，用异常捕获处理中间可能出现的截断。核心逻辑不复杂，但有一个点很容易被忽略：真实雷达文件偶尔会有脏字节或者长度标记缺失，这个时候如果直接崩掉，整个批处理就停了。我选择用while加异常捕获的循环结构，读到异常就跳过这一条记录并记录日志，而不是直接退出。

def parse_volume(filepath): with open(filepath, 'rb') as f: volume = RadarVolume() while True: try: rec_len = struct.unpack('<I', f.read(4))[0] except struct.error: break if rec_len == 0: break body = f.read(rec_len - 4) if len(body) < rec_len - 4: break tail = struct.unpack('<I', f.read(4))[0] if tail != rec_len: continue parse_record(body, volume) return volume

这个函数在真实项目中跑过，稳定性尚可。有一个处理技巧值得分享：用while加异常遍历记录结构，比一次性读全文件再切分更稳健，因为脏字节经常出现在你意想不到的位置。parse_record函数内部根据记录类型再次分发，目前我支持了反射率、径向速度、速度谱宽三类基础产品，以后增加差分相移或者相关系数时，只在分发处加一个分支即可，不用动主体逻辑。

4.3 性能优化与批量处理

当文件很多时，解码性能会成为瓶颈。一个台站一天数百个时次，每个时次几十兆字节，如果解码逻辑写得太慢，整个分析流程就会被拖垮。优化第一思路是尽量用numpy整列操作，避免逐库for循环。比如把一条径向的原始数据一次性读成uint16数组，然后整体转成float32并做换算，比每条数据循环处理快一个量级。

第二条是使用numba对纯数值循环做JIT加速。不过在我的场景里，解码部分已经足够快，主要耗时反而在matplotlib绘图上。实际经验是，把单个时次的解码时间从最初的几十秒优化到几秒，主要靠两条：将径向数据先读成uint16数组再整体换算，以及解码结果缓存到内存，批处理画图时同一个体扫对象不解码两次。批量处理还可以引入multiprocessing并行，按文件号拆分任务，多核并行解码。如果只是画图，不建议过度设计，顺序处理在十几分钟级别的任务量上完全够用。

5. 常见问题与排查技巧实录

5.1 解码阶段常见报错

让我列一张排查速查表，这些都是我在实际调试中遇到的问题，不一定每次都会遇到，但遇到之后能快速定位方向。

解码阶段遇到乱码，先别急着改代码。我习惯用十六进制查看器打开文件，对比已知格式文档里的固定字节，确认文件头里的魔数和站点编号是否对上。这个习惯能省下大量排查时间，尤其是在你刚拿到一个新地区或者新型号雷达数据的时候。文件头匹配通过后，再继续往下读径向数据，定位问题就快很多。

5.2 绘图阶段常见问题

绘图阶段的问题通常是视觉上的，比如画面出现明显的空洞、颜色不对、图太大跑不动等。画面出现空洞，多为无效值nan在插值时产生的空白区域，可以在插值前把nan替换成背景值或使用nearest插值。颜色与业务不一致，检查色标边界值是否设置正确，比如反射率图上出现负回波渐变，说明vmin设得不对。出图文件过大，多半是网格分辨率设得太高，PPI网格分辨率设成500乘以500已经足够，再高意义不大；如果确实需要高分辨率，优先提高数据端分辨率而不是绘图网格。

还有一个很容易被忽略的点是坐标轴刻度。雷达图通常需要保持纵横比相等，也就是set_aspect('equal')，否则圆形探测范围会被拉伸成椭圆，方位角视觉上就会失真。叠加地图时也要注意投影统一，matplotlib的Basemap虽然老旧，但cartopy现在用得更普遍，如果不想引入额外依赖，只画简单的经纬网和边界线也可以。

5.3 环境配置建议

环境方面，我建议Python 3.8以上，依赖包就四个：numpy、matplotlib、scipy、pytest，其中pytest是可选的。安装命令很简单：

pip install numpy matplotlib scipy pytest

如果遇到matplotlib中文字体不显示，需要在plot_radar.py开头配置中文字体。这个问题在科研绘图中非常常见，标题里的中文和负号容易变成方框或者异常显示。

import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS'] plt.rcParams['axes.unicode_minus'] = False

做了这个配置后，标题、图注里的中文才能正常显示，负号也不会变成竖条。建议在项目一开头就配好，不要等到画了几十张图之后才发现中文全变方框。

我在实际调试这套代码时最大的体会是：雷达数据处理的难点不在写代码，而在先吃透格式文档。CINRAD的二进制格式里埋着不少历史兼容字段，不同厂家、不同型号会有细微差异，拿着官方说明一字节一字节对照着读，比随便搜一段代码来改要靠谱得多。绘图部分则多试几个配色，出了第一张能看的图之后，后面就是批量出图的体力活了。

这个工具后续还可以扩展的方向很多，比如接入批量下载任务、把解码结果输出成NetCDF或GTiff、增加切片和剖面交互式查看。建议先把底层解析和绘图流程跑通，再做增量功能，每一步都能实实在在往上叠，这也是我维护这个项目一路过来的方法。

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