☰
PyCINRAD实战:从天气雷达基数据到PPI图像绘制全流程
2026/10/2 2:46:04 网站建设 项目流程

一年多以前我第一次拿到一份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本身,而是绘图侧的cartopy

PyCINRAD 安装本身没什么难度:

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 python=3.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, projection=ccrs.PlateCarree()) mesh = ax.pcolormesh(lon, lat, data, cmap=cmap, norm=norm, shading='auto', transform=ccrs.PlateCarree()) ax.coastlines(linewidth=0.8) ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False) cb = plt.colorbar(mesh, ax=ax, extend='both', shrink=0.8) cb.set_label('Reflectivity (dBZ)') plt.title('PPI 0.5 deg') plt.savefig('ppi_example.png', dpi=150, bbox_inches='tight')

这里有几个细节值得展开。

雷达站位置我建议单独标注。上面代码里,雷达站经纬度可以从数据文件属性里取,不同版本属性名不一样,有的叫经纬度数组,有的还要从文件头解析。最简单的方法是在图上画一个三角形标记代表雷达站位置,直接从lon[0,0]和lat[0,0]取。严格来说lon[0,0]是第一个径向库第一个距离库的经纬度,通常就是雷达站附近,用于标注完全够用。

shading='auto'这个参数也很关键。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里的transform=ccrs.PlateCarree()表示数据经纬度本身是WGS84坐标,但axes的投影也是一个等经纬度投影,两层叠起来看似一致,可如果海岸线数据源本身带精度偏移,或者画图范围跨的纬度比较大,PlateCarree的线性经纬度网格会造成距离变形。

我的处理思路是:业务展示图若无特殊要求,固定用PlateCarree+range设置,保证出图可复现;如果做中纬度强对流分析,推荐将axes投影换成LambertConformal,强回波区域的形状更接近实际:

proj = ccrs.LambertConformal(central_longitude=120.0, central_latitude=30.0, standard_parallels=(30.0, 60.0)) ax = fig.add_subplot(111, projection=proj)

这里有一个容易忽略的点:改axes投影后,ax.pcolormesh里的transform参数仍然要写ccrs.PlateCarree(),因为雷达数据本身是经纬度坐标,这一步是告诉cartopy数据点的原始坐标系,和axes投影不是一回事。

4.2 单位坑:dBZ和线性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 色标选择:为什么不能随手用rainbow

matplotlib自带的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, cmap=cmap, norm=norm, shading='auto') plt.colorbar(mesh, ax=ax, shrink=0.8) plt.title(f.split('_')[4][:12] + ' UTC') plt.savefig(f'frame_{i:03d}.png', dpi=100) images.append(Image.open(f'frame_{i:03d}.png')) plt.close() images[0].save('radar_animation.gif', save_all=True, append_images=images[1:], duration=300, loop=0)

有几个细节我踩过。第一,循环里必须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层仰角)时,用缓存方案从原来的十几分钟压缩到了两分钟以内。这个优化思路,在数据量上来之后非常值得做。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询