雷达基数据这玩意儿第一次拿到手的人大概率会愣在原地说不出话一个几十MB的二进制文件没有图像预览没有标准抬头你甚至不知道它属于哪个雷达站、几点几分扫的。我这些年用Python做气象雷达数据处理从最早的逐字节解码到现在固定依赖PyCINRAD最大的感受就是“把基数据变成PPI图”这件事真正难的不是画图那一下而是前面的格式解码、坐标变换和数据清洗。这篇文章就把我跑通全流程的路径完整记录下来从雷达基数据长什么样开始讲到用PyCINRAD读出数据、理解数据结构最后输出一张能直接用于分析的PPI图。适合两类人看一是刚接触雷达数据、连基数据和产品图都分不清的新手二是已经装了PyCINRAD但卡在“读不出来”或者“出图不对”的开发者。1. 先搞清基数据里到底装了什么再动手写代码1.1 体扫、仰角、距离库这三个概念决定你怎么读数据雷达基数据本质上是一次“体积扫描”的原始记录。天气雷达的天线在一个固定仰角下360度旋转一边转一边发射电磁脉冲并接收回波这个动作叫“体扫”的一部分。一次完整体扫通常包含多个仰角从低到高依次覆盖例如0.5°、1.5°、2.4°……具体层数由雷达的扫描模式VCP决定业务上常见的是9层或14层左右。每一层里雷达在某个方位角上对指定距离范围内连续采样采样点的编号叫“距离库”不同的最大探测距离对应不同的距离库长度比如230公里范围通常对应几百到上千个库。当你打开一个基数据文件里面实际存的就是这些层、这些方位角、这些距离库上的原始测量值。不同物理量会分别存储反射率单位是dBZ速度单位是m/s谱宽单位也是m/s。对于双偏振雷达还有差分反射率、差分相移、相关系数等衍生观测。换句话说文件的结构是一个“仰角层 × 方位角 × 距离库 × 物理量”的四维档案PyCINRAD给我们的Radar对象不过是把这个四维档案按层、按物理量切出来后的二维视图。1.2 文件名里的信息别浪费拿到文件第一步先别急着往代码里塞路径花十秒钟看下文件名。S波段新型基数据常见命名方式是这样的Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin拆开看Z9898是雷达站站号20230723090000是观测时间年月日时分秒依次排列SA表示雷达型号是SA型DOR表示这是一部多普勒雷达CAP代表数据格式版本。有些老数据是.RD后缀命名风格类似站号和时间位置都在文件名里。这些信息非常有用尤其是在批量处理时不需要打开文件就能根据文件名筛选时次、筛选站点。更重要的是遇到解析失败时文件名直接提示了雷达型号你可以有针对性地去查这个型号的格式手册而不是拿着文件一头雾水。1.3 为什么基数据不能直接当图像处理这里补充一个最容易被忽略的认知问题。很多教程会让你把雷达数据“读进数组然后画图”但如果你直接对原始二进制做这个操作大概率得到一堆没有意义的数字。原因有三层。第一二进制编码问题。雷达文件的头结构、字段长度、字节序在不同厂家里是不完全一样的必须按对应雷达型号的格式规范去解析否则读出来全是“错的数字”。第二数据网格问题。雷达记录的是极坐标网格方位角和距离才是它的天然坐标轴直接当普通图像矩阵处理几何关系从第一步就错了。第三质量控制问题。原始基数据里存在大量缺测值、超出量程的异常值、地物杂波、速度模糊等瑕疵这些如果不做标记或过滤画出来的图会把噪声当回波。PyCINRAD存在的意义就是把这些最繁琐又最基础的工作统一封装好让我们能把精力放在上层分析上。2. Python侧环境准备装好PyCINRAD这一步问题最多2.1 安装依赖链与最常翻车的地方PyCINRAD的安装命令很简单一行pip install cinrad但这一行背后有几个依赖容易让人崩溃。它依赖numpy、scipy、matplotlib、cartopy、netCDF4、pyproj、cftime等一堆库其中cartopy属于比较重的依赖在Windows上直接用pip装可能会因为缺少GEOS等底层库而报错。我的建议是如果你电脑上还没有科学计算环境先装一个Anaconda或Miniconda然后用conda把cartopy装好再回来pip install cinrad。如果网络慢把PyPI源或者conda源换成国内镜像能省大量时间。已经装了Python但还没装过科学计算库的直接pip install numpy scipy matplotlib cartopy cinrad会把整条依赖链一起拉下来但注意看安装日志里有没有编译错误。注意安装完成后第一件事是在Python里执行import cinrad如果这一步报错大概率是某个依赖库缺失或版本冲突先把报错发到搜索引擎解决完再继续。不要一上来就处理数据环境不干净会污染后面所有排查。2.2 没有雷达数据怎么练手用官方样例数据没有雷达数据是练不了手的PyCINRAD在io模块里提供了一个下载样例数据的功能可以直接把网上的标准格式基数据拉回本地。实际用法很简单一个函数调用传入一个文件名它会自动下载到当前目录。我第一次用的时候就想早知道有这个功能前面就不用为到处找数据发愁了。如果你是刚接触这一套强烈建议先下载一两个不同型号的样例文件把读取、画图的流程全部跑通再去折腾自己单位或公开渠道拿到的正式数据。还有一个好处是样例数据的时间、站号都是固定的出问题可以比对官方文档的示例结果排查起来比用自己的数据方便得多。import cinrad filename cinrad.io.load_standard_data( Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin ) print(filename)2.3 装完先验货跑出第一张图再想别的环境装好、样例数据下载好之后先别急着深入细节直接跑一个最简脚本确认整条链路是通的。这段代码如果环境正确几十秒内就会生成一张反射率PPI图import cinrad from cinrad.io import CinradReader from cinrad.visualize import PPI filename cinrad.io.load_standard_data( Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin ) f CinradReader(filename) r f.get_data(1, 230, REF) # 第1层仰角230公里范围反射率 fig PPI(r, cmapref) fig(first_ppi.png)这里的get_data三个参数各有用意第1个参数是仰角索引一般从1开始对应体扫里的第一层第2个是最大距离范围单位公里第3个是物理量类型REF代表反射率。PPI类负责整套坐标变换和底图渲染fig(first_ppi.png)是保存图片的入口。跑通这一张图后面所有问题都值得继续往下看。3. 读取基数据两条读取路径和一个关键数据结构3.1 新旧格式的读取入口怎么选PyCINRAD根据基数据文件的格式主要提供两个读取入口。新格式文件文件名带CAP、SAD等标识一般是.bin或.dat用CinradReader老格式文件常见.RD后缀用StandardData。在新版本里CinradReader的兼容性已经做得不错很多老格式也能解析但你拿到的文件一旦出现读取异常第一反应应该是检查你选没选对入口而不是怀疑库有问题。两个类的用法基本一致先实例化传入文件路径然后调用get_data方法拿到某个仰角、某个物理量的数据对象。如果你不知道自己文件属于哪种可以先打开文件头部几个字节看看有没有可识别的ASCII标识或者干脆两个类各试一次能正常读出来并返回非空数据就算过。3.2 理解返回对象的关键属性get_data返回的对象关键属性就那几个搞明白之后你就有了一副看数据的“眼镜”。属性含义常用说明data物理量数据二维numpy数组shape为(径向数, 距离库数)azimuth方位角数组每个径向对应一个方位角单位度range距离库数组每个距离库的斜距单位米el当前仰角值单位度方便你确认读的是哪层site站点信息包含雷达站号、经纬度等scan_time体扫开始时间批量归档时非常有用dtype物理量类型字符串比如REF、VEL这里最关键的是要明白data的坐标系是极坐标横坐标是距离库索引纵坐标是方位角索引直接print出来的数组看起来会像一张扇形的“原始图”但这还远不是一张合格的PPI图。很多卡住的同学就是栽在这一步数据读出来了也print出来了但不知道下一步该怎么把它变成正常图像。3.3 拿到数据先自检三个指标读出来的数据不要直接画图先做三个小检查。第一用np.isnan看看无效值占比如果数据全是NaN多半是仰角参数或物理量类型写错了或者文件本身缺层。第二看data.min()和data.max()反射率的合理范围一般在-30到80dBZ之间如果出现极端异常值比如上千说明数据解析可能有字节序问题或者这个文件根本不是基数据。第三检查azimuth数组的长度正常应该和data.shape[0]一致如果不一致说明文件可能损坏。这三个检查两分钟内做完能避免后面排查两个小时的尴尬。3.4 一个稳妥的读取函数示例把上面串起来一个稳妥的读取函数可以写成这样import numpy as np from cinrad.io import CinradReader, StandardData def read_radar(filepath, tilt1, drange230, dtypeREF): # 根据文件内容选择读取入口 with open(filepath, rb) as fh: head fh.read(8) readers [CinradReader, StandardData] for Reader in readers: try: f Reader(filepath) r f.get_data(tilt, drange, dtype) if r is not None and not np.all(np.isnan(r.data)): return r except Exception: continue raise ValueError(file not parseable: {}.format(filepath))这个函数属于个人习惯不是标准API但它很稳定地帮我自动判断新老格式。注意with open只读了个头部就关闭了不会占用文件句柄太久。4. PPI图背后的坐标算法从极坐标到平面坐标4.1 一层的基数据就是一个极坐标网格PPI图的全称是Plan Position Indicator中文常叫平面位置显示图。它的意思是把雷达在某个固定仰角扫描得到的回波投影到一个水平面上用颜色表示回波强度。雷达在扫描时记录的数据是按“方位角”和“距离”两个维度组织的这就是极坐标网格。以一层的反射率数据为例每一条径向对应一个方位角从0度到360度每条径向上又有一串距离库代表从雷达到远处的一段段距离。整个二维数组的行是方位角列是距离。要把它显示成一个正常的平面图最直观的做法是把每个格点从极坐标(r, θ)换算成平面坐标(x, y)再把每个格点的物理量值填到对应位置。4.2 坐标映射其实就是高中那点三角函数极坐标转平面的公式大家高中学过x r·sin(θ)y r·cos(θ)。雷达方位角的约定是正北方向为0度顺时针增加所以正东是90度正南是180度。这里r要用实际水平距离θ要用弧度制。低仰角扫描时仰角本身带来的斜距和水平距离差异很小工程上经常直接用斜距作水平距离处理但如果你做的是精确计算或者仰角超过10度就必须乘cos(仰角)把斜距投影到水平面。具体搬进数组里你需要在每个(distance, azimuth)对上做换算。假设距离库是range单位米径向的方位角是azimuth单位度那么对于一个径向i它的X坐标是range * sin(radians(azimuth[i]))Y坐标是range * cos(radians(azimuth[i]))。把所有格点都这样换算一遍就得到一张平面散点分布每个点带一个回波值再经过网格插值把空白区域填满这就是PPI图的核心算法。4.3 插值与扇形盲区为什么图上会有空缺极坐标点转换到直角坐标后分布并不是均匀填满整个方框的。靠近雷达的地方点数密集远处点数稀薄而且如果雷达的扫描模式存在某个仰角的遮挡还会出现整片扇形盲区。所以直接scatter画点会显得很挤、很多空洞。PyCINRAD在这块内部做了插值处理把极坐标数据重采样到规则的矩形网格上遇到缺测或超出最大探测距离的位置留空。理解这一点对排查图像问题很重要你看到PPI图上有个扇区是空白的不一定是你代码有bug可能是雷达本身的扫描盲区也可能是距离库截断之后插值区域的边界效应。4.4 一个常被忽略的细节PPI不是严格的水平面还有一个容易被忽略的专业细节天气雷达的波束是倾斜向上扫描的加上大气折射作用波束路径会缓慢弯曲。同一张PPI图上靠近雷达的位置回波高度比较低到边缘位置波束高度可能已经到几公里以上。也就是说PPI图严格来说是“准水平面”产品它把不同高度的回波投影到一个平面上了。日常看强对流回波时大家默认这个误差可以接受但你如果要叠加地面观测、做定量降水估计就需要考虑到边缘位置的回波并不代表地面附近的真实情况。这个认知能帮你避免很多“明明PPI图上回波很强地面却没下雨”的困惑。5. 实战从基数据到一张能直接用的PPI图5.1 最简版本四行代码出图最简版本四行就能出图from cinrad.io import CinradReader from cinrad.visualize import PPI f CinradReader(Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin) r f.get_data(1, 230, REF) # 第1层仰角230km距离库反射率 PPI(r)(output.png)如果你是从头开始建议先跑通这一版确认数据、环境都没有问题再往里面加东西。PPI类在实例化时就会根据数据的dtype决定默认色标和量程PPI(r)(output.png)这个写法看着有点绕其实第一条语句创建了图对象第二条语句就是把图像保存到文件。保存成功后你会看到当前目录多了一张PNG图上是一个以雷达为中心的圆形区域颜色代表反射率。5.2 定制版换物理量、调色标、控制范围默认图可用但很多时候你想调整色标、配色和范围。PPI类支持的常见参数包括cmap、norm等。cmap可以直接传类似ref、vel、sw这样的物理量色标名也可以连matplotlib里通用的颜色映射。实际项目中我习惯显式传cmap避免不同版本默认值变化导致出图风格漂移。范围的调整一般通过get_data里的drange参数控制比如我只关心150公里以内的回波就把drange改成150图中心区域的细节会被放大边缘也被裁掉。这里要注意drange的最大值受文件本身距离库长度的限制你填500但如果文件本身就只有230公里的探测范围实际还是以文件为准。from cinrad.io import CinradReader from cinrad.visualize import PPI f CinradReader(Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin) for tilt in [1, 2, 3]: r f.get_data(tilt, 150, REF) PPI(r, cmapref)(ppi_tilt{}.png.format(tilt))5.3 叠加站点和地名沿用matplotlib的思路单看色斑图不够做业务报告时一般需要叠加雷达站点位置、周边城市名、省界这些要素。PyCINRAD的PPI对象在创建时会尝试用Cartopy绘制底图所以你先出来的一张图里通常已经有海岸线和省界了这一点比很多从零手写的方案省事。如果需要再叠加城市名等文本标注可以在PPI对象的基础上继续操作它内部保存的matplotlib坐标轴。不过不同版本的内部属性名略有差异建议打开IDE补全看下你装的那个版本暴露了哪些属性。如果你对底图的样式要求特别高我的经验是干脆不用PyCINRAD的默认底图只用它把色斑图画好然后导出数据再用cartopy自己重新绘制整张图可控性最强。5.4 批量出图一个体扫生成多张PPI业务上很少只画一张图更多是一次处理一个体扫的所有仰角、或者一个时段里多个体扫的固定仰角。批量处理要注意两点。第一get_data每调用一次都会重新解析一次数据如果你要对同一个文件取多个物理量、多个仰角尽量在循环外层保留读取对象不要每层都重新new一个CinradReader。第二如果批量处理大量文件建议给文件名加上站点、时间、仰角、物理量这些关键标识统一输出规范后面做时间序列分析时能省很多事。from pathlib import Path from cinrad.io import CinradReader from cinrad.visualize import PPI input_file Path(Z_RADR_I_Z9898_20230723090000_O_DOR_SA_CAP.bin) out_dir Path(ppi_output) out_dir.mkdir(exist_okTrue) f CinradReader(str(input_file)) for tilt in range(1, 10): try: r f.get_data(tilt, 230, REF) except Exception: print(tilt {} read failed.format(tilt)) continue png out_dir / {}_tilt{:02d}_REF.png.format(input_file.stem, tilt) PPI(r, cmapref)(str(png))这个循环会对一个体扫里前9层仰角各出一张反射率PPI图如果某层本身不存在捕获异常跳过而不是让整个批量任务中断。6. 我踩过的坑和对应的排查思路6.1 整张图一个颜色先别怀疑软件最常遇到的问题就是图出来了但整张图几乎一个颜色没有回波结构。先别怀疑软件按顺序排查第一确认你选的物理量名是对的把REF写成REFL这种错误很常见get_data内部一般会抛异常或返回异常结果第二检查返回的r.data是不是全为NaN如果是说明该层仰角在该距离范围内没有有效回波这种时候图自然是一片空白第三检查drange和文件实际探测范围是不是差太远你把230公里的数据硬要画成500公里边缘大部分都是无效值图看起来也会很空第四如果以上都没问题看看是不是cmap选了个不合适的配色导致强回波和弱回波颜色区分度很低。6.2 图上回波位置“歪”了多半是坐标系问题图上回波位置和真实地理方位对不上也是常见的坑。先检查雷达坐标来源PyCINRAD读取时用的是文件头里的站点经纬度如果文件头本身经纬度异常整张图对地图的位置就全偏了。第二确认你的底图是同一个坐标系。PyCINRAD默认的底图经过投影转换你如果自己往上面叠加普通经纬度点坐标系不统一就会偏。第三有时候看起来“歪”其实是图的纵横比问题PNG图的宽高比例和经纬度跨度不一致导致的视觉变形尤其当你在小范围切图时会出现把图像尺寸设置成和范围比例一致就能缓解。6.3 中文地名全部变成方框中文字体问题叠加中文标注时最常见的坑是matplotlib默认字体不支持中文出图后名字全部变成小方框。解决方法也不复杂在绘图前把matplotlib默认字体切换成系统里支持中文的字体比如SimHei、Noto Sans CJK等这一步在Windows、Linux、macOS上略有差异但思路都一样先把字体注册给matplotlib再设置rcParams。如果你在服务器上跑服务器通常没有中文字体更省事的方案是把标注改成拼音或者英文避免字体问题拖慢出图进度。6.4 批量处理时内存和速度先想清楚再循环基数据文件动辄几十上百MB批量处理时内存很容易爆。经验上有两个实用技巧。一是适度缩小drange230公里和460公里距离库数差一倍插值计算量会差好几倍二是不要一次性把所有仰角、所有物理量都读进内存再统一画图而是“读一层、画一层、释放一层”。如果你的分析只需要反射率一个物理量就只调REF不要好奇地把VEL、SW都拉进来数组加对象本身都很占内存。批量几百个文件时建议用multiprocessing做并行但并行进程数不要超过CPU核心数太多否则IO反而会成为瓶颈。6.5 出图清晰度不够保存参数里藏玄机最后补一个不算坑但很影响体验的细节出图清晰度。默认保存的PNG分辨率对屏幕看完全足够但如果要放进报告里打印建议在保存时提高dpi。PyCINRAD的PPI对象生成的是一个matplotlib Figure你可以用自带保存接口也可以在拿到Figure对象后调用它的savefig传dpi参数。更进阶的做法是把回波数据先导出成标准格式比如NetCDF再用统一的绘图脚本生成整套报告图这样出图风格完全可控也和PyCINRAD的版本更新解耦。我自己的经验是雷达基数据处理这件事工具箱里PyCINRAD只占了很小一部分真正值钱的是你对自己数据的了解。你花一下午把某个站点的文件结构、常见坏数据特征、扫描模式都摸清楚之后所有画图、分析都是水到渠成。最后分享一个自动化小技巧如果你每天都要定时处理新到的基数据可以把流程封装成一个脚本用系统自带的任务计划程序或cron定时执行生成的PPI图加上站号、时次、仰角的规范命名直接进业务网页展示这套东西一旦跑起来能省下大量重复手工劳动。