简介基于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或FortranPython开发周期短很多对科研场景尤其友好而且兼容性广泛从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, errorsignore)这里需要注意格式化串里的小于号表示小端千万别漏掉。之前有同事拿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, dtypenp.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), methodlinear) zi np.ma.masked_invalid(zi) fig, ax plt.subplots(figsize(8, 6)) pm ax.pcolormesh(xi, yi, zi, cmapreflectivity, vmin0, vmax70) plt.colorbar(pm, axax, labeldBZ) 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剖面与多仰角拼接除了PPIRHI在分析回波垂直结构时非常常用。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 解码阶段常见报错让我列一张排查速查表这些都是我在实际调试中遇到的问题不一定每次都会遇到但遇到之后能快速定位方向。问题现象可能原因解决方案struct.unpack报参数长度不足偏移量或记录长度算错检查记录长度标记用记录头偏移调试读出的站点号是乱码字节序或偏移错误确认小端字节序对照格式文档偏移表距离库数量为0或异常大读取到非径向数据块增加记录类型判断过滤无效记录换算后的物理量全是极值无效值标记未处理设置合理的无效值上下限替换为nan解码阶段遇到乱码先别急着改代码。我习惯用十六进制查看器打开文件对比已知格式文档里的固定字节确认文件头里的魔数和站点编号是否对上。这个习惯能省下大量排查时间尤其是在你刚拿到一个新地区或者新型号雷达数据的时候。文件头匹配通过后再继续往下读径向数据定位问题就快很多。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、增加切片和剖面交互式查看。建议先把底层解析和绘图流程跑通再做增量功能每一步都能实实在在往上叠这也是我维护这个项目一路过来的方法。本文还有配套的精品资源点击获取