简介面向气象数据处理与雷达技术学习者的多普勒雷达强度分析工具包围绕敏视达雷达回波数据展开可直接运行查看降水强度分布。程序基于多普勒频移原理完成距离解析、速度计算和强度评估配套真实探测数据适合大气科学专业学生、气象工程师及雷达算法初学者用于理解强度处理流程并练习基础图像显示。压缩包共13个文件、约737KB核心包含可执行exe及对应cpp源程序便于对照研读算法逻辑数据文件以dat为主另含调试与工程配置类文件整体结构简洁。已有257人学习下载借助该资源可快速掌握雷达强度数据显示方法为后续分析反射率因子、判别强对流天气打下实操基础。1. 从 .rar 里的“intensity”说起气象雷达强度图到底能干什么拿到一份标注为 ShowRadarData_intensity.rar 的气象雷达资料包绝大多数人的第一反应是解压、找脚本、跑一下然后对着黑黢黢的控制台发愣。这个包的核心价值不在压缩包本身而在“intensity”这个词——它是多普勒气象雷达最基本的观测产品反射率强度场。无论你手里的雷达是 SA 波段还是 CB 波段无论数据格式是国产雷达基数据还是通用二进制产品最后都要落到一件事把一串二进制数字变成一张能看出回波强弱的平面图。这篇文章只讲一件事——怎么从“拿到包”走到“画出图”中间的数据格式、坐标换算、色标、投影和各类翻车现场,一次说透。2. 多普勒雷达强度数据的真实结构为什么 ShowRadarData 要先过格式与坐标两关很多人把“气象雷达数据处理”想成读文件、画图两步实际跑起来才发现中间横着两座山一是数据格式千差万别二是极坐标数据不能直接往直角坐标的画布上摆。ShowRadarData 这类小工具诞生就是为了翻过这两座山。2.1 强度场不是一张照片距离库、方位角与扫描策略多普勒天气雷达在体扫模式常见的是 VCP21 这类扫描策略下天线按多个仰角做 360° 旋转。每转一圈雷达在一个仰角上按固定方位角间隔比如 1°采样每个方位角上沿径向按固定距离间隔记录能量回波强度。这样得到的强度场天然是一张极坐标“扇形图”——横轴是方位角纵轴是距离数据点均匀分布在极坐标网格上而不是你想当然的横平竖直的像素网格。这里的“距离库”指每个距离采样点它的间隔由脉冲宽度决定常见有 250m、1km 两档。以 250m 库长、230km 最大探测距离为例一根径向只有 920 个距离库但整个体扫会包含几十个仰角、每个仰角 360 根径向——数据总量并不小。把这样一组数据还原成图像第一件要做的事是理解每个数值在物理空间里的位置它属于哪个仰角、哪个方位角、哪个距离库。这三个坐标缺一个画出来的图就是废的。强度场的物理含义是反射率因子单位是 dBZ。不同强度的回波对应不同降水过程10dBZ 以下的通常是层状云弱回波35dBZ 以上基本可以判定为中等以上对流降水50dBZ 以上则可能伴随冰雹。这就是为什么很多脚本里会硬编码一组阈值——不是随便拍的是气象业务里约定俗成的分级。提示读数据前先确认雷达的扫描模式。有的数据是单仰角 PPI有的是一次完整体扫的多个仰角。ShowRadarData 这类工具通常默认输出某个指定仰角的 PPI 图不要在体扫数据上只取一层就算完事。2.2 先探格式再动手识别文件头、仰角表与每一个“径向”我拿到一份未知格式的雷达基数据第一件事不是写完整解析器而是用十六进制编辑器或 Python 直接看文件头。常见的气象雷达基数据文件会以一定长度的文件头开头里面写着站点名称、雷达型号、体扫开始时间、仰角层数、每层的仰角值和径向数。文件头之后是一个接一个的径向数据块每个径向数据块通常包含方位角信息、该径向的质量控制标记、以及一串强度数据。用一个通用的读取框架做“格式侦察”比一上来就套某个标准库要稳妥。下面这段代码作用是打开一个未知二进制文件观察文件头区域的前若干字节找出可打印字符——它们往往就是站点名和时间字符串。import struct def probe_radar_header(filepath, max_bytes1024): 探测雷达基数据文件头找出可打印字符和疑似时间字段。 适合在拿到未知 .bin / .dat 文件时做第一轮判断。 with open(filepath, rb) as f: raw f.read(max_bytes) # 把前 1024 字节按一行 16 字节排版打印出来 for offset in range(0, len(raw), 16): chunk raw[offset:offset16] hex_part .join(f{b:02X} for b in chunk) ascii_part .join(chr(b) if 32 b 127 else . for b in chunk) print(f{offset:04X} {hex_part:48} {ascii_part})这段探测代码逻辑很简单每次取 16 字节左边打印十六进制字节右边打印可打印字符的 ASCII 形式。跑完之后你会看到类似“Z9090”“2024-07-05”“VCP21”这样的字段从乱码中浮出来。站点编号、时间、扫描模式一旦定位再去查对应雷达的数据格式说明就有方向了。从文件头拿到关键元数据之后第二步是定位径向数据块的起始偏移。常见做法是固定文件头长度或者从文件头解析出“数据体起始字节数”。我习惯把文件头解析结果直接打印出来然后把指针移动到数据区起点逐根径向读取方位角、仰角和强度长度字段。下面这段代码演示了最通用的逐径向读取逻辑它假设每根径向有一个 2 字节的方位角单位度可能带小数、一个 2 字节的径向长度表示该径向有多少距离库之后紧跟着等长的强度数据。import numpy as np def read_radial_intensity(filepath, header_len256, max_range_bins1000): 按“方位角(2字节) 库数(2字节) 强度数据(N字节)”的通用结构 读取一个仰角文件中的所有径向返回方位角列表和强度矩阵。 如果实际格式不同通常会在这里暴露字节错位导致强度值全是乱数。 with open(filepath, rb) as f: f.seek(header_len) # 跳过文件头 raw f.read() # 剩余部分全部读入 azimuths [] intensities [] offset 0 while offset 4 len(raw): az struct.unpack_from(H, raw, offset)[0] / 100.0 # 低字节在前按 1/100 度存储 offset 2 n struct.unpack_from(H, raw, offset)[0] offset 2 if n max_range_bins or n 0: break # 库数异常说明字节错位或已读到尾部 data np.frombuffer(raw, dtypenp.uint8, countn, offsetoffset) offset n azimuths.append(az) intensities.append(data) return np.array(azimuths), np.array(intensities)这段代码的关键参数有三个header_len是文件头长度max_range_bins是距离库数量上限用来挡住字节错位时的异常值方位角按“低字节在前”的小端序读取这是 Windows 系统生成的数据最常见的形式。如果你的文件是 Sun 服务器上转出来的老数据可能要用H读大端序跑出来方位角全在 0 到 655 之间乱跳第一步就应该怀疑字节序。注意很多基数据格式里坏库值不是 0 而是极值比如 65535 或者把强度存成 0~255 的索引值而不是直接存 dBZ。解析时先看强度数组的最大值和最小值再决定要不要做查表映射。常见做法是文件里存的只是强度等级索引真正的 dBZ 值要查一个 256 级线性表。2.3 仰角、时次、雷达型号ShowRadarData 这类工具“吃”什么元数据解析器真正依赖的元数据没几项仰角层的总数、当前层在体扫里的序号、站点经纬度画底图要用、雷达最大探测距离和库长。很多脚本把“站点经纬度”写死在某张配置文件里换站点就出错——这是最常见的低级翻车现场。正确做法是先从文件头读读不到再落到配置。我一般会在解析时把这四样东西一并打印出来体扫开始时间、当前仰角、径向数量、最大强度值。这四个值对不上后面画图基本是白画。3. 把强度场从极坐标搬到平面核心换算与最小可跑脚本格式解析过了接下来是把极坐标的强度数据映射到水平的笛卡尔网格上。这个过程叫坐标换算市面上有 pyart、wradlib 这类现成库做这件事但自己写一遍的意义在于你能在出错时知道是哪一环的问题。3.1 距离库索引怎么变成平面坐标斜距、水平距离与地球曲率修正雷达波束不是贴着地面走的。从天线出发的波束按仰角向上倾斜加上地球是圆的波束的“高度”随距离增长——这就是经典的雷达方程几何关系。把一个距离库索引换算成水平面上的(X, Y)至少要经过三步索引乘库长得斜距 R用仰角算出该斜距对应的水平投影距离 S再修正地球曲率带来的高度误差。如果不做曲率修正在 230km 处的高度偏差可能达到几公里回波位置会明显偏移。最常见的简化公式是用“4/3 地球半径”模型把地球半径等效放大三分之一用于估算波束中心高度而水平距离直接用斜距乘以仰角余弦来近似。对于 20° 以下的小仰角这种近似在工程上完全够用。换算输出的是一个二维网格上每个像素的经纬度或本地笛卡尔坐标。3.2 一个能跑通的最小脚本读完一个仰角并映射到二维网格下面这段代码把上一节读出的强度矩阵映射到笛卡尔网格上。核心思路先创建目标平面网格网格点间距设为 1km覆盖范围由最大显示距离决定然后对每个网格点算出对应的方位角和距离库索引从极坐标强度矩阵里取值。为了避免空值点采用最近邻取值不插值——因为气象回波本身就存在大量无效区强行平滑会制造假回波。import numpy as np def polar_to_cartesian(azimuths, intensity_matrix, range_res250, max_range230000, grid_size230): 将极坐标强度数据映射到 (grid_size x grid_size) 的笛卡尔网格。 azimuths: 每个径向的方位角单位度 intensity_matrix: 形状 (num_radials, num_bins) 的强度值 range_res: 距离库长米 max_range: 最大显示距离米 返回水平和垂直坐标网格、强度网格。 az_rad np.deg2rad(azimuths) # 方位角转弧度 num_radials, num_bins intensity_matrix.shape # 建立目标平面-max_range 到 max_range步长为 range_res axis np.arange(-max_range, max_range range_res, range_res) X, Y np.meshgrid(axis, axis) dist np.sqrt(X**2 Y**2) # 每个网格点到原点的距离 angle np.arctan2(Y, X) * 180.0 / np.pi # 网格点方位角单位度 # 把方位角归一到 0-360 angle (angle 360.0) % 360.0 # 距离索引网格距离 / 库长越界标记为 -1 dist_idx np.floor(dist / range_res).astype(int) dist_idx[dist_idx num_bins] -1 # 方位角归一到 0-360并找最近的径向索引 ang_idx np.floor((angle - azimuths[0]) / (360.0 / num_radials)).astype(int) ang_idx (ang_idx num_radials) % num_radials grid np.full_like(X, np.nan, dtypefloat) valid dist_idx 0 grid[valid] intensity_matrix[ang_idx[valid], dist_idx[valid]] return X, Y, grid这段代码里的参数有讲究angle np.arctan2(Y, X)得到的是逆时针方向角如果雷达数据的方位角是“北零度、顺时针”两者需要再校准一次方向关系ang_idx做了取模运算保证方位角 359° 和 0° 能正确衔接dist_idx对超远距离做了静默截断避免索引越界报错。最后的掩膜方式是这套流程决定成败的一步把所有无效距离库先压在 NaN 上再交给绘图库处理不要用 0 填充——否则一张回波图会平白多出一大圈“零回波”假象。注意这段代码默认雷达数据是由正北方向开始、按顺时针方位角排列。国产 S 波段雷达和进口雷达在这点上常常相反最直观的验证方法是找一张已知回波形态的个例看回波移动方向是否合理。3.3 径向缺测与坏库的掩膜时机插值前必须处理的脏数据雷达在工作时会有径向数据缺失某根方位角上没有任何返回、某几根径向数据被遮挡地物杂波抑制过度这些脏数据如果直接参与最近邻取值会在图上形成放射状亮线或空洞。我的做法是在映射循环之前先做一次“径向健康度检查”统计每根径向的有效库数有效库数低于整根径向 10% 的直接把该径向强度全部置为无效。坏库值通常在整个数据文件里是统一的比如 -999、65535、255 这类标记。读取完成后立刻统一掩膜不要等画图时再处理。掩膜后的强度矩阵配合上一节的最近邻映射生成的网格上无效区域呈现自然的空洞后续无论是叠加地图还是做动画都不会出现恼人的白色噪点。4. 画一幅能看出门道的反射率图色标、投影与产品输出坐标换算完成数据已经在直角网格上了接下来的问题是“怎么画才不出错”。色标选错会误导人投影选错图就废了。4.1 反射率色标为什么不能随手选台风、层状云与强对流的 dBz 分级气象业务里反射率色标有一套约定俗成的分段方式10dBZ 以下用冷色深蓝、青10-35 用绿到黄35-50 用橙到红50 以上用紫甚至洋红。这套配色不是审美问题是有可读性要求的——业务人员只需扫一眼颜色分布就能判断出哪里有强对流中心。反过来说如果随手用 matplotlib 默认的 viridis 或者 jet回波强度的梯度感会被完全带偏弱回波和强回波之间的过渡在视觉上会被夸大或抹平。所以我在自己的脚本里一直维护着一张写死的颜色查找表对应 -10 到 70dBZ每 5dBZ 一档。外面的库画图再方便我也只在“快速预览”时用默认色标出正式产品一定切回这张表。4.2 叠加地理信息的标准做法用 Cartopy 给强度图配底图单纯画强度场看不出回波和地形、城市、河流的关系做灾害天气复盘时必须叠加地理信息。Cartopy 是目前最稳的选择它能直接读自然地球要素或者本地 shp 文件。下面这段代码展示如何把强度网格和海岸线叠加import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature def plot_ppi(X, Y, grid, lon_ref, lat_ref, output_path): 在指定经纬度中心点绘制 PPI 强度图叠加海岸线和省界。 X, Y: 平面上距离雷达的笛卡尔坐标网格单位米 grid: 强度网格单位 dBZ无效区为 NaN lon_ref, lat_ref: 雷达站点经纬度 crs ccrs.PlateCarree() fig plt.figure(figsize(10, 10)) ax fig.add_subplot(1, 1, 1, projectioncrs) # 把以雷达为原点的平面坐标转成经纬度 # 简化的等距换算一度纬度约 111 km lon lon_ref X / (111000 * np.cos(np.deg2rad(lat_ref))) lat lat_ref Y / 111000.0 cmap, norm build_dBZ_colormap() # 业务色标自动归一化 img ax.pcolormesh(lon, lat, grid, cmapcmap, normnorm, transformcrs) ax.add_feature(cfeature.COASTLINE, linewidth0.5) ax.add_feature(cfeature.BORDERS, linewidth0.5) ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()]) plt.colorbar(img, shrink0.8, labeldBZ) plt.savefig(output_path, dpi150, bbox_inchestight)这段代码里最容易出错的是经纬度换算把笛卡尔距离转经纬度时经度的单位长度会随纬度变化不能用一个固定系数。111000 * np.cos(np.deg2rad(lat_ref))是常用的近似做法用于局部小范围时误差可忽略如果整个图覆盖范围超过 200km建议用更严格的投影方式如兰勃特正形投影。build_dBZ_colormap是我自己封装的色标生成函数它会根据数据范围自动选择从哪一档开始填色避免一张弱回波图被默认色标渲染成暴雨级别。4.3 输出 PNG 还是 GeoTIFF按报告、网页与业务存档三种场景选产品图输出格式取决于用途。出报告或放在网页上PNG 足够注意 dpi 不能低于 150否则印刷会糊做空间分析或业务存档一般输出 GeoTIFF把经纬度信息直接写进文件头后续可以用 GIS 软件直接叠加载入。还有一份常用输出是每帧带时间戳的序列图文件名里必须包含年月日时分否则做动画回放时排序就乱了。下表是我平时在不同场景下的输出参数用途格式分辨率/参数备注汇报材料PNGdpi200RGB叠加站点名和时次不要叠加行政边界网页展示PNGdpi100压缩级别 6控制体积黑白底图可以压缩业务存档GeoTIFF原分辨率无压缩必须写入雷达站点经纬度和投影信息数据分析NetCDF浮点保留 NaN后续用 Python/R 读取最方便输出文件名命名习惯建议与此保持一致{site}_{yyyymmddHHMM}_{elev}.png这样排序时天然按时间排列。5. 四个避坑记录与解析失败排查从 NaN 到“雷达图长反了”读雷达数据经常会翻车以下四条是实践中最常见的高频问题每一条都按“现象 → 原因 → 解决”的思路来说。5.1 解压即崩编码问题之外还要小心“rar 用来加载广告的子程序”现象从网上下载的 ShowRadarData_intensity.rar 解压后双击脚本报错显示找不到某个文件或乱码。更糟的一种情况是解压过程中弹出了不明窗口或浏览器广告脚本根本没法正常跑。原因这类资料包经过多次转手文件名可能带着 GBK 编码的中文目录Windows 解压工具在默认编码下会把路径搞乱还有一种情况是包里被塞进了“rar用来加载广告的子程序”这样的捆绑程序解压工具一并释放了出来。解决用 Python 的 zipfile 或 7-Zip 指定 UTF-8 解码重新打包解压解压后先看包里是否有可疑 exe看到直接删然后检查默认编码是否为 UTF-8再跑脚本。资料包从公共渠道下载后尤其是经过了几道转手的运行前一定要检查文件哈希或勒索病毒特征不能盲目双击。5.2 全图 NaN坏库标记没掩膜仰角索引还读错现象脚本跑通了输出的强度网格全是 NaN图上一片全空。原因雷达数据里的坏库值并不是 0而是 65535 或者 255 这类哨兵值如果没有做掩膜映射时会把哨兵值当作有效强度算进去然会在数值上出现极端值整个色标的 colorbar 被拉伸到失真另一种可能是读仰角时索引取错读出来的是噪声层或空层。解决读取强度矩阵后立刻做一次极值检查打印最大值和最小值如果是哨兵值用np.where把它们统一替换成 NaN然后把每个仰角的“有效库数占比”打印出来占比低于 5% 的层直接丢弃不要参与绘图。5.3 回波左右颠倒方位角零点的起点和扫描方向没有对齐现象强度图画出来了但回波的形态是对称反转的就像照了镜子或者整块回波被旋转了一个固定角度。原因雷达数据的方位角有两种排列习惯——正北方起始、顺时针扫描行业常用和正东方起始、逆时针扫描少数进口雷达。如果映射脚本里假设了其中一种遇到另一种就必然翻车。解决先取任意一个包含明显回波特征的个例打开一张已知正确的回波图做对比然后在映射前检查方位角数组的第一个值和最后一个值如果方位角从 359° 递减到 0°应按相反方向重新索引再不行就把网格上的 X、Y 互换或者翻转排掉左右镜像问题。5.4 老资料包忘记解压密码先找注释再别传播现象下载的课程资料类 rar 包或者同事转发的老雷达数据包解压时提示输入密码密码没人记得。原因打包者当年设了密码防篡改时间久了密码失传这是“课程资料.rar 忘记解压密码”这类问题的高发场景。解决先用7z l -slt查看压缩包注释很多打包者会把密码写在注释里再试试文件名、站点号、日期等常规信息组合如果都失败不要花大量时间去暴破除非你确认这些数据对该项目不可替代——常见做法是重新找数据源或者找原始设备的记录。这类包解不开时最忌讳的是到处转手这不仅浪费别人的时间还可能传播被破解过的恶意程序。提示拿到任何老雷达数据包建议第一时间用 7-Zip 重新打包一次并去掉旧密码或加密头避免后续接手的人再次被挡在门外。5.5 程序在 32 位环境上读取超过 2GB 数据时内存爆炸现象脚本在处理长时间序列或全仰角体扫数据时内存占用飙到几个 GB然后程序被系统杀掉。原因一次性把所有仰角的强度矩阵读进内存再同时保留极坐标和笛卡尔两份数据内存轻松翻倍。解决改为逐仰角读取、映射、释放只保留当前结果不保留中间数组如果实际项目需要同时处理几十个时次就把输出直接写成 NetCDF 而不是留在内存里。内存换时间这条在雷达数据处理上基本不划算——一份体扫数据动辄上百 MB做动画回放时逐帧处理才是正路。6. 把强度帧串成动画做回放验证解析结果的一个可重复技巧当你第一次把静态图跑通时不要马上宣布大功告成。真正的验证方法是把这个站多个时次的强度图串成动画观察回波的移动是否连续、是否有跳变、旋转方向是否正确。这一步能同时暴露 5.3 节的镜像问题和 3.3 节的径向缺测问题如果回波移动方向和实际盛行风方向相反或者回波像“跳帧”一样突然消失再出现说明前面的某一步还有隐患。动画生成不需要额外装太多依赖Matplotlib 的FuncAnimation配合PIL保存 PNG 序列然后合成 GIF 是最省事的做法。我一般把每帧间隔设为 200ms对应 6 分钟一个体扫的数据看起来节奏刚好。import os import numpy as np from PIL import Image def make_animation_from_frames(frame_paths, output_gif, duration_ms200): 把按时间排序的 PNG 帧序列合成为 GIF。 frame_paths: 按时间顺序排列的文件路径列表 duration_ms: 每帧播放时长200ms 对应 6 分钟一个体扫的常见节奏 frames [Image.open(p) for p in frame_paths] # 统一尺寸避免帧与帧之间大小不一致导致 GIF 抖动 sizes set(f.size for f in frames) if len(sizes) 1: target max(sizes) frames [f.resize(target) for f in frames] frames[0].save(output_gif, save_allTrue, append_imagesframes[1:], durationduration_ms, loop0, optimizeFalse)这段代码省掉了读数据的环节直接把画好的静态图合成动画。参数上只有duration_ms需要按数据时间间隔调整5 分钟一个体扫就设 250ms10 分钟一个体扫就设 400ms。optimizeFalse是有意为之——气象图细节多开启优化会压缩掉弱回波区域的颜色渐变不利于肉眼识别。我自己的习惯是任何一次新数据源接入都先用一段 2 小时的连续帧动画跑一遍验收再把静帧放入业务汇报材料。这一关过了后面的数据分析、算法评估才真正站得住脚。希望这套流程能帮你少走一段弯路。本文还有配套的精品资源点击获取