简介:面向气象数据处理与雷达技术学习者的多普勒雷达强度分析工具包,围绕敏视达雷达回波数据展开,可直接运行查看降水强度分布。程序基于多普勒频移原理完成距离解析、速度计算和强度评估,配套真实探测数据,适合大气科学专业学生、气象工程师及雷达算法初学者用于理解强度处理流程并练习基础图像显示。压缩包共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_bytes=1024): """ 探测雷达基数据文件头:找出可打印字符和疑似时间字段。 适合在拿到未知 .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:offset+16] 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_len=256, max_range_bins=1000): """ 按“方位角(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, dtype=np.uint8, count=n, offset=offset) 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_res=250, max_range=230000, grid_size=230): """ 将极坐标强度数据映射到 (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, dtype=float) 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, projection=crs) # 把以雷达为原点的平面坐标转成经纬度 # 简化的等距换算:一度纬度约 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, cmap=cmap, norm=norm, transform=crs) ax.add_feature(cfeature.COASTLINE, linewidth=0.5) ax.add_feature(cfeature.BORDERS, linewidth=0.5) ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()]) plt.colorbar(img, shrink=0.8, label='dBZ') plt.savefig(output_path, dpi=150, bbox_inches='tight')这段代码里最容易出错的是经纬度换算:把笛卡尔距离转经纬度时,经度的单位长度会随纬度变化,不能用一个固定系数。111000 * np.cos(np.deg2rad(lat_ref))是常用的近似做法,用于局部小范围时误差可忽略;如果整个图覆盖范围超过 200km,建议用更严格的投影方式(如兰勃特正形投影)。build_dBZ_colormap是我自己封装的色标生成函数,它会根据数据范围自动选择从哪一档开始填色,避免一张弱回波图被默认色标渲染成暴雨级别。
4.3 输出 PNG 还是 GeoTIFF:按报告、网页与业务存档三种场景选
产品图输出格式取决于用途。出报告或放在网页上,PNG 足够,注意 dpi 不能低于 150,否则印刷会糊;做空间分析或业务存档,一般输出 GeoTIFF,把经纬度信息直接写进文件头,后续可以用 GIS 软件直接叠加载入。还有一份常用输出是每帧带时间戳的序列图,文件名里必须包含年月日时分,否则做动画回放时排序就乱了。下表是我平时在不同场景下的输出参数:
| 用途 | 格式 | 分辨率/参数 | 备注 |
|---|---|---|---|
| 汇报材料 | PNG | dpi=200,RGB | 叠加站点名和时次,不要叠加行政边界 |
| 网页展示 | PNG | dpi=100,压缩级别 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_ms=200): """ 把按时间排序的 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_all=True, append_images=frames[1:], duration=duration_ms, loop=0, optimize=False)这段代码省掉了读数据的环节,直接把画好的静态图合成动画。参数上只有duration_ms需要按数据时间间隔调整:5 分钟一个体扫就设 250ms,10 分钟一个体扫就设 400ms。optimize=False是有意为之——气象图细节多,开启优化会压缩掉弱回波区域的颜色渐变,不利于肉眼识别。
我自己的习惯是,任何一次新数据源接入,都先用一段 2 小时的连续帧动画跑一遍验收,再把静帧放入业务汇报材料。这一关过了,后面的数据分析、算法评估才真正站得住脚。希望这套流程能帮你少走一段弯路。
本文还有配套的精品资源,点击获取