一年多以前我第一次拿到一份CINRAD/SA天气雷达基数据时对着满屏二进制字节完全没辙反射率、速度、谱宽几大字段挤在同一个文件里手动解析到凌晨才勉强读出一段数据。后来同事甩给我一个库名——PyCINRAD我才发现原来绘制雷达PPI图像这件事本来就是几行代码的工夫。这篇博客就把我实际跑通的经验和踩过的坑完整写出来特别是坐标投影、色标、单位换算这类的隐蔽问题。内容适合气象专业学生、做短临预报分析的研究人员以及想用Python处理雷达数据的入门爱好者。1. 为什么最终选了PyCINRAD雷达基数据这块硬骨头的正确打开方式1.1 自己解析雷达基数据会遇到的三堵墙国内新一代多普勒天气雷达的基数据文件名常常长这样Z_RADR_I_Z9898_20230704120000_O_DOR_SA_CAP.bin。这类数据本身是二进制格式SA/SAD、CB、SC等雷达型号之间还有差异不同型号的字节对齐方式、径向库数量、分层策略完全不同。直接解析的话第一堵墙就是格式文档难找且版本混乱第二堵墙是径向数据的极坐标投影到经纬度网格需要自己处理地球曲率第三堵墙是就算画出来了色标、距离环、雷达站标注这些细节也能磨掉一整晚。我最早尝试过自己写解析脚本确实能出图但代码又臭又长换一个批次的文件经常要重新调参。后来也试过一些重量级气象库但为了画一张PPI去搭整套依赖性价比太低。PyCINRAD恰好卡在功能足够用和上手足够轻的位置上。1.2 PyCINRAD和同类库的取舍下面这个对比表是我实际折腾过之后整理的直接说结论方案上手难度对国内雷达数据兼容性可定制化程度适合场景纯手写解析高只适配单一型号完全可控研究算法底层Py-ART中高需要转换器格式映射偏弱强科研级质控与算法开发wradlib中高提供CINRAD reader但依赖较重强学术研究与教学PyCINRAD低原生支持SA/SAD等国内常见格式中等快速可视化、业务成图、入门教学我的选择逻辑很直接先跑通业务再考虑精调。PyCINRAD 把读取、投影、绘图封装得足够好所谓的投影背后也已经实现了标准大气条件下的球面坐标换算不需要我手撸四参数公式。等到后面需要做定量分析和质量控制时再回头对接Py-ART也不迟。2. 环境准备与数据读取检查画图之前先避开两个暗坑2.1 安装依赖重点不是PyCINRAD本身而是绘图侧的cartopyPyCINRAD 安装本身没什么难度pip install pycinrad如果网络不好用国内镜像pip install pycinrad -i https://pypi.tuna.tsinghua.edu.cn/simple真正的坑在绘图侧。PyCINRAD 提供了快速绘图接口但如果你希望出图时叠加海岸线、行政边界这些地图要素绕不开 cartopy。cartopy 的底层依赖 GEOS、PROJ 这些C库直接用 pip 装经常出现版本不匹配表现是装完导入时报一个 DLL load 相关的错或者画地图时投影报错。我的建议是直接用 conda 创建独立环境conda create -n radar python3.10 conda activate radar conda install cartopy pip install pycinrad这样cartopy的C库依赖交给conda处理能省掉很多莫名其妙的问题。Python版本建议3.9到3.11之间太新的Python版本有时会让某些科学计算包还没跟上。2.2 文件名、路径和编码看着不起眼实际最常卡住雷达基数据文件名里通常包含雷达站号、观测时间、雷达型号等信息本身是ASCII字符。但如果你是找气象台拷贝的数据拷出来经常变成中文名或者带空格的名字比如7月4日雷达数据.bin。PyCINRAD读取中文路径有时会碰编码问题建议先把数据统一改成英文名放到纯英文路径下再处理mv 7月4日雷达数据.bin Z_RADR_I_Z9898_20230704120000_O_DOR_SA_CAP.bin数据放好后读取和基本检查也很简单import pycinrad radar pycinrad.io.radar(Z_RADR_I_Z9898_20230704120000_O_DOR_SA_CAP.bin) print(radar)不同小版本的API可能略有差异有的版本写的是pycinrad.io.radar()新一些的版本可能是pycinrad.io.read_cinrad()本质一样。装好后我不建议死记API直接在Jupyter里对radar对象按Tab补全或者用dir(radar)看一下可用的属性和方法比查文档更快。2.3 体扫数据合法性抽查先看仰角序列再决定画哪层雷达基数据是体扫一个文件里有多层仰角的扫描数据。拿到文件后第一件事不是急着画图而是先打印扫描信息确认这个体扫完不完整、各仰角层数据有没有缺# 查看扫描信息不同版本属性名可能有差异 try: print(radar.scan_info) except AttributeError: for i, angle in enumerate(radar.angles): print(i, angle)打印出来你会看到类似0.5度、1.5度、2.4度、3.4度...这样的仰角序列。这一步非常关键因为如果体扫只进行到一半就中断后面层是全空的画出来就是一张缺了半边数据的图。3. PPI图像绘制核心代码从基数据对象到一张业务级可用的图3.1 快速出图先看效果再谈优化PyCINRAD 自带可视化接口最快的方式是直接用它的封装方法import pycinrad import matplotlib.pyplot as plt radar pycinrad.io.radar(Z_RADR_I_Z9898_20230704120000_O_DOR_SA_CAP.bin) # 画第0层仰角的PPI fig, ax plt.subplots(figsize(10, 10)) pycinrad.visualize.plot_ppi(radar, 0) plt.show()这一版图能出来但说实话距离业务级还有距离。自带的封装在色标选择上通常比较随意图上也没有地图要素对外的正式图件基本不会直接用。但这一步用来快速检查数据质量、确认回波范围效率很高。我的习惯是先跑这版确认数据没问题再进手动绘制流程。3.2 手动绘制PPI可控的色标、坐标和地图叠加业务成图或者写论文插图建议手动控制每个环节流程拆开是四步取网格数据、过滤噪声、定义色标、叠加地图要素。import pycinrad import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from matplotlib.colors import BoundaryNorm, ListedColormap radar pycinrad.io.radar(Z_RADR_I_Z9898_20230704120000_O_DOR_SA_CAP.bin) # 取第0层仰角、230km范围、反射率产品 ppi radar.get_ppi(0, 230, REF) lon ppi[lon] lat ppi[lat] data ppi[data] # 过滤掉弱回波噪声业务上经常把低于0 dBZ的部分置为无效 data np.ma.masked_where(data 0, data) # 线性Z转dBZ的兜底判断防止单位坑 if np.nanmax(data) 100: data 10 * np.log10(data)色标采用气象上常用的NWS风格从灰色到绿色再到黄橙红分段之间用BoundaryNorm切分levels [-20, -10, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70] colors [#000000, #7e7e7e, #00b4e6, #01d300, #00e600, #2ab800, #73ce00, #c8fe00, #fefe00, #ffce00, #ffa800, #fe8400, #ff6800, #ff0000, #c80000, #a00000, #7a0000] cmap ListedColormap(colors) norm BoundaryNorm(levels, cmap.N) fig plt.figure(figsize(10, 10)) ax fig.add_subplot(111, projectionccrs.PlateCarree()) mesh ax.pcolormesh(lon, lat, data, cmapcmap, normnorm, shadingauto, transformccrs.PlateCarree()) ax.coastlines(linewidth0.8) ax.gridlines(draw_labelsTrue, dmsTrue, x_inlineFalse, y_inlineFalse) cb plt.colorbar(mesh, axax, extendboth, shrink0.8) cb.set_label(Reflectivity (dBZ)) plt.title(PPI 0.5 deg) plt.savefig(ppi_example.png, dpi150, bbox_inchestight)这里有几个细节值得展开。雷达站位置我建议单独标注。上面代码里雷达站经纬度可以从数据文件属性里取不同版本属性名不一样有的叫经纬度数组有的还要从文件头解析。最简单的方法是在图上画一个三角形标记代表雷达站位置直接从lon[0,0]和lat[0,0]取。严格来说lon[0,0]是第一个径向库第一个距离库的经纬度通常就是雷达站附近用于标注完全够用。shadingauto这个参数也很关键。get_ppi返回的经纬度网格和数值网格有时候维度不匹配如果不加这个参数pcolormesh会报维度不一致的错。自动shading会处理好边界坐标这个细节。3.3 版本差异引出的API使用经验我这边用的版本radar.get_ppi(0, 230, REF)三个参数分别表示第几层仰角、最大距离km、产品类型。但如果你安装的版本API有变化不要慌帮助文档和源码都能救你help(radar.get_ppi)PyCINRAD有一段时间API调整比较频繁比如get_ppi在新版本里的函数签名可能换成了别的方式返回的字典字段也可能从data变成了别的名字。我的经验是看返回字典的keys()ppi radar.get_ppi(0, 230, REF) print(ppi.keys())如果你的版本返回的字段跟我这里不一样按实际字段名调整就行。核心思想是雷达基数据最终要落到三个量——经度网格、纬度网格、反射率数值网格拿到这三个量画图就成功了一大半。4. 避坑指南坐标误差、色标陷阱与杂波干扰的完整排查链路4.1 坐标投影为什么叠加地图后回波位置会偏有一次我在图上加了海岸线发现强回波的中心位置和自动站雨量对不上偏差大概有十几公里。第一反应以为是雷达标定问题查了半天发现是坐标投影的理解出了偏差。PPI本质上不是平面扫描而是雷达波束在固定仰角上绕垂直轴旋转形成的锥面扫描。把锥面上的数据点投影到平面地图时必须处理两个问题一是波束斜距到地面水平距离的换算二是地球曲率和大气折射引起的波束高度抬升。PyCINRAD内部默认采用标准大气折射模型也就是等效地球半径取实际半径的4/3倍这是大多数天气雷达业务软件的标准假设。那偏差是怎么来的我遇到的情况是叠加地图时投影坐标系没有对齐。pcolormesh里的transformccrs.PlateCarree()表示数据经纬度本身是WGS84坐标但axes的投影也是一个等经纬度投影两层叠起来看似一致可如果海岸线数据源本身带精度偏移或者画图范围跨的纬度比较大PlateCarree的线性经纬度网格会造成距离变形。我的处理思路是业务展示图若无特殊要求固定用PlateCarreerange设置保证出图可复现如果做中纬度强对流分析推荐将axes投影换成LambertConformal强回波区域的形状更接近实际proj ccrs.LambertConformal(central_longitude120.0, central_latitude30.0, standard_parallels(30.0, 60.0)) ax fig.add_subplot(111, projectionproj)这里有一个容易忽略的点改axes投影后ax.pcolormesh里的transform参数仍然要写ccrs.PlateCarree()因为雷达数据本身是经纬度坐标这一步是告诉cartopy数据点的原始坐标系和axes投影不是一回事。4.2 单位坑dZ和线性Z的分不清问题雷达产品里的反射率因子有两种表达方式线性值Z单位mm^6/m^3和分贝值dBZ。两者的关系是dBZ 10 * log10(Z)我拿到过一批数据画出来的图整体颜色特别奇怪回波区全部是深红色色标的数值区间也不对从几百到几万。后来排查发现get_ppi返回的这个字段是线性Z值不是dBZ。业务上大家默认看的dBZ所以必须转换。上面代码里if np.nanmax(data) 100这个判断就是用来兜底的因为正常dBZ数值区间在-20到70之间超过100的时候说明数据可能是线性Z值。更好一点的做法是直接看ppi.keys()里有没有单位字段或者打印数据的统计量print(np.nanmin(data), np.nanmax(data))如果最小值是负数那基本可以确定是dBZ因为线性Z不可能是负值。这也是快速判断单位的一个技巧。4.3 地物杂波和超折射回波一张干净图背后的数据质量控制天气雷达最低仰角0.5度的PPI在晴空条件下经常能看到雷达站附近的零散杂波尤其是雨天超折射现象会把地物回波放大表现为从站址向外辐射的条状、瓣状回波。这些杂波在图上会误导判断。我常用的简单处理方式是# 距离库太近的区域容易受地物杂波影响但不直接置零只做展示层过滤 data np.ma.masked_where(data 0, data)把0 dBZ以下的背景噪声过滤掉。这个操作只适合出展示图千万不能认为这就是质量控制后的定量数据。如果要正儿八经做定量降水估算需要做多普勒速度退模糊、杂波抑制、衰减订正等一系列处理建议到时改用Py-ART这类库来做。还有一个小技巧对比不同仰角的PPI。地物杂波通常只出现在最低仰角且近距离范围内如果1.5度仰角同一位置没有回波最低仰角那里的强回波大概率是杂波可以辅助人工判断。4.4 色标选择为什么不能随手用rainbowmatplotlib自带的jet或者viridis画天气雷达图颜色和回波强度没有统一语言。写过气象报告的人都知道业务上要求同一种回波强度在所有图里都显示成同一种颜色这样横向对比多时次图像才可靠。NWS风格色标是这个领域的默认实践从灰到绿再到黄橙红对应回波从弱到强。如果你不想自己定义颜色列表也可以直接尝试导入一些现成库import pyart cmap pyart.graph.cm_colorblind_reflectivity不过为了不引入额外重量级依赖我自己更习惯把颜色表写成一个独立的py文件需要时直接import。上面代码里的颜色列表就是我从NWS标准色标整理出来的覆盖了-20到70 dBZ的常规区间。5. 进阶玩法批量成图、地图投影选择与动画输出5.1 多时次批量出图一次跑完一个降水过程业务分析很少只看单张PPI我经常要处理连续两三个小时的体扫数据每次体扫一个文件批量出图再拼动画能直观看到回波的生消演变import glob import numpy as np import matplotlib.pyplot as plt from PIL import Image files sorted(glob.glob(data/20230704*.bin)) images [] for i, f in enumerate(files): radar pycinrad.io.radar(f) ppi radar.get_ppi(0, 230, REF) data np.ma.masked_where(ppi[data] 0, ppi[data]) fig, ax plt.subplots(figsize(8, 8)) mesh ax.pcolormesh(ppi[lon], ppi[lat], data, cmapcmap, normnorm, shadingauto) plt.colorbar(mesh, axax, shrink0.8) plt.title(f.split(_)[4][:12] UTC) plt.savefig(fframe_{i:03d}.png, dpi100) images.append(Image.open(fframe_{i:03d}.png)) plt.close() images[0].save(radar_animation.gif, save_allTrue, append_imagesimages[1:], duration300, loop0)有几个细节我踩过。第一循环里必须plt.close()不关闭的话内存会被不断增长的Figure对象吃掉第二文件名的排序要用sorted(glob.glob(...))否则frame_010会排在frame_001前面造成时间轴错乱第三如果要处理几百个文件建议改成多进程批量绘图但多进程里不要直接调用plt.show()只做保存。5.2 投影不是越复杂越好LambertConformal确实在中纬度看起来舒服但有一个副作用在跨度过大的图上lat/lon网格线会弯得厉害色块边缘也容易变形。做单站雷达图时我通常先在PlateCarree下看整体形态需要出正式图再切Lambert。实际上对单张PPI来说PlateCarree在200km半径范围内产生的形变肉眼很难察觉不需要为了显得专业而强行上复杂投影。5.3 PPI的局限和CAPPI扩展思路最后聊一个容易忽略的物理本质PPI是固定仰角的锥面扫描波束高度随距离增大而抬升。0.5度仰角在100km处波束中心高度大约在1.5到2公里这就意味着远处看到的回波并不是地面附近的降水而是中低层的降水结构。所以分析层状云降水时不能只看最低仰角PPI就下结论。如果要更准确地反映某一高度层的降水分布应该用CAPPI等高平面位置显示它的思路是把多个仰角的数据插值到同一个固定高度上。PyCINRAD对CAPPI也有部分支持但插值算法相对简单结果精度和Py-ART或wradlib比要粗糙一些。我的建议是日常快看用PyCINRAD出PPI深入研究降水结构时换专业库做CAPPI。实测下来PyCINRAD最适合它的场景就是快速制图和业务初判。一张PPI图从读文件到保存最多也就几秒时间中间最耗时间的反而是等get_ppi做坐标变换。文件多的时候建议把转换好的经纬度网格先存成npy缓存第二次读取直接加载能省掉一大半时间。我处理一次完整降水过程48个时次、每时次4层仰角时用缓存方案从原来的十几分钟压缩到了两分钟以内。这个优化思路在数据量上来之后非常值得做。