文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >专栏 >雷达系列 | Py-ART 多雷达组网实战（grid_from_radars）

雷达系列 | Py-ART 多雷达组网实战（grid_from_radars）

作者头像
用户11172986
发布于 2026-05-13 17:26:49
发布于 2026-05-13 17:26:49
950
举报
文章被收录于专栏：气python风雨气python风雨

雷达系列 | Py-ART 多雷达组网实战（grid_from_radars）

上周朋友丢过来三部菲律宾雷达的 CAPPI 文件，问能不能拼成一张完整的回波图。单部雷达半径只有几百公里，覆盖不全。这种事 Py-ART 一行命令就能搞定，但参数要怎么配、网格范围设多大才不至于一片空白，背后还是有些细节值得说一说。

本篇是雷达系列教程的 Py-ART 分篇，专注 pyart.map.grid_from_radars() 这一条链路。wradlib 那条更灵活的链路放在另一份单独的教程里，两份内容相互独立，可分开学习。

本教材基于 Py-ART 2.2.0 版本，早期版本中 grid_from_radars 的参数命名（如 roi_func 与 weighting_function）存在差异，遇到报错先确认版本号。

阅读路径

  1. 数据介绍
  2. 读取多部雷达
  3. GateFilter 数据质控
  4. grid_from_radars 网格化组网
  5. Cartopy 出版级绘图
  6. 完整代码与附录
代码语言：javascript
复制
!pip install arm-pyart --upgrade -i https://pypi.mirrors.ustc.edu.cn/simple/

1. 数据介绍

本教程使用 3 部菲律宾地区的雷达 CAPPI 数据，文件格式为 CF/Radial NetCDF。

雷达站

文件名

纬度

经度

数据形状 (azimuth, range)

SBAG (Baguio)

20251103051105_SBAG_CAPPI_2km.nc

16.356°N

120.559°E

(360, 480)

SBAL (Baler)

20251103051001_SBAL_CAPPI_2km.nc

15.749°N

121.632°E

(360, 480)

SDAET (Daet)

20251103051038_SDAET_CAPPI_2km.nc

14.129°N

122.983°E

(360, 960)

CAPPI（Constant Altitude Plan Position Indicator）即等高度位显示产品，本数据为 2 km 高度的 CAPPI 反射率因子。

数据中的反射率场名称为 Reflectivity_CAPPI_2km，单位 dBZ。三部雷达分布在吕宋岛北、中、南三处，覆盖范围互有重叠，正好适合演示组网。

代码语言：javascript
复制
# 导入常用库
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import pyart
import cmaps
import warnings
warnings.filterwarnings('ignore')

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

file_baguio = '20251103051105_SBAG_CAPPI_2km.nc'
file_baler  = '20251103051001_SBAL_CAPPI_2km.nc'
file_daet   = '20251103051038_SDAET_CAPPI_2km.nc'

print('依赖检测完成')

2. Py-ART 组网

Py-ART（Python ARM Radar Toolkit）是 ARM 计划维护的雷达数据处理工具包，其中 pyart.map.grid_from_radars() 函数把多部雷达数据直接插值到统一的笛卡尔网格上，一步完成组网。整体流程拆成读数据、做质控、跑网格化、出图四个动作。

2.1 读取多部雷达数据

用 pyart.io.read_cfradial 直接读 CF/Radial 格式的 NetCDF。读完后得到三个 pyart.core.Radar 对象，里面是按极坐标存储的扫描数据。

代码语言：javascript
复制
# 使用 pyart.io.read_cfradial 读取 CF/Radial 格式的 NetCDF 文件
radar_baguio = pyart.io.read_cfradial(file_baguio)
radar_baler  = pyart.io.read_cfradial(file_baler)
radar_daet   = pyart.io.read_cfradial(file_daet)

# 查看雷达基本信息
for name, radar in [('SBAG', radar_baguio), ('SBAL', radar_baler), ('SDAET', radar_daet)]:
    print(f'雷达站， {name}')
    print(f'  字段， {list(radar.fields.keys())}')
    print(f'  位置， ({radar.latitude["data"][0]:.3f} N, {radar.longitude["data"][0]:.3f} E)')
    print(f'  数据形状， {radar.fields["Reflectivity_CAPPI_2km"]["data"].shape}')
    print()

2.2 数据质控（GateFilter）

组网前先做一遍质控，剔除两类东西。一类是扫描过渡区（transition）也就是天线在仰角切换瞬间采到的奇怪数据，另一类是明显异常的超高值，比如某些原始文件里会出现 Reflectivity > 200 dBZ 的填充值。Py-ART 的 GateFilter 把这两类直接打成 mask，后续插值时自动忽略。

实际工作里要不要应用质控，看个人偏好。原始数据本身已经比较干净时，可以先不传 gatefilter 看一眼组网效果，再决定要不要加。

代码语言：javascript
复制
def make_gatefilter(radar):
    """创建门控过滤器，排除过渡扫描和异常高值。"""
    gf = pyart.filters.GateFilter(radar)
    gf.exclude_transition()                              # 排除扫描过渡区
    gf.exclude_above('Reflectivity_CAPPI_2km', 100)      # 排除超过 100 dBZ 的异常值
    return gf

# 分别为三部雷达创建过滤器
gatefilter_baguio = make_gatefilter(radar_baguio)
gatefilter_baler  = make_gatefilter(radar_baler)
gatefilter_daet   = make_gatefilter(radar_daet)

print('GateFilter 创建完成')

2.3 网格化与组网（grid_from_radars）

pyart.map.grid_from_radars() 是 Py-ART 组网的核心函数，主要参数列在下面这张表里。

参数

说明

grid_shape

网格形状 (z, y, x)，本例只组网 2 km 高度，z 维度取 1

grid_limits

网格范围 ((z_min, z_max), (y_min, y_max), (x_min, x_max))，单位米，以雷达阵列中心为原点

fields

需要插值的字段列表

weighting_function

权重函数，'barnes2' 表示 Barnes 二阶权重

gridding_algo

网格化算法，'map_gates_to_grid' 把雷达门直接映射到网格

roi_func

影响半径函数，'dist_beam' 根据波束距离自动计算

gatefilters

门控过滤器元组（可选），用于数据质控

注意，grid_shape 中的 y/x 维度越大，分辨率越高，但计算时间也越长。做演示用 300×300 即可，做正式产品再上 800×800。

代码语言：javascript
复制
# 执行网格化组网
grid = pyart.map.grid_from_radars(
    (radar_daet, radar_baguio, radar_baler),                # 多部雷达元组
    # gatefilters=(gatefilter_daet, gatefilter_baguio, gatefilter_baler),  # 可选，应用质控
    grid_shape=(1, 300, 300),                                # (z, y, x) 网格尺寸
    grid_limits=((0, 2000),                                  # 垂直方向 0~2000 m
                 (-2000000.0, 2000000.0),                    # 南北方向 ±2000 km
                 (-1000000.0, 1000000.0)),                   # 东西方向 ±1000 km
    fields=['Reflectivity_CAPPI_2km'],                      # 反射率场
    weighting_function='barnes2',                            # Barnes2 权重函数
    gridding_algo='map_gates_to_grid',                       # 网格映射算法
    roi_func='dist_beam',                                    # 距离波束影响半径
)

print('组网完成')
print(f'网格字段， {list(grid.fields.keys())}')
print(f'网格形状， {grid.fields["Reflectivity_CAPPI_2km"]["data"].shape}')

2.4 提取并绘制组网结果（Cartopy）

组网完成后，从 Grid 对象里取出第 0 层（也就是 2 km CAPPI 这一层）的反射率，和对应的经纬度二维数组。Py-ART 的 Grid 已经内置了 point_latitude / point_longitude，不需要自己做坐标转换。

代码语言：javascript
复制
# 从 Grid 对象中提取经纬度和反射率数据
refl = grid.fields['Reflectivity_CAPPI_2km']['data'][0]   # 取第 0 层（2 km 高度）
lats = grid.point_latitude['data'][0]
lons = grid.point_longitude['data'][0]

print(f'反射率范围， {np.nanmin(refl):.2f} ~ {np.nanmax(refl):.2f} dBZ')
print(f'纬度范围， {lats.min():.2f} ~ {lats.max():.2f}')
print(f'经度范围， {lons.min():.2f} ~ {lons.max():.2f}')
代码语言：javascript
复制
# 使用 Cartopy 绘制组网反射率图（出版级质量）
cmap = cmaps.radar_1
level = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]

# 设置出版级字体
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman', 'DejaVu Serif']
plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams['font.size'] = 11
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['axes.titlesize'] = 13
plt.rcParams['xtick.labelsize'] = 10
plt.rcParams['ytick.labelsize'] = 10

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 8), dpi=300)
ax = plt.axes(projection=proj)
ax.set_extent([120, 126, 12, 16.5], crs=proj)

# 绘制反射率填色
mesh = ax.contourf(
    lons, lats, refl,
    transform=proj,
    cmap=cmap,
    levels=level,
)

# 海岸线
ax.coastlines(resolution='10m', linewidth=0.6, color='dimgray')

# 网格线，淡化处理
gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.3,
    color='gray',
    alpha=0.4,
    linestyle='--',
    xlocs=range(116, 129, 2),
    ylocs=range(12, 21, 2),
)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}

# 雷达站位置，空心圆 + 白色描边，避免遮挡数据
radar_coords = {
    'SBAG':  (float(radar_baguio.latitude['data']), float(radar_baguio.longitude['data'])),
    'SBAL':  (float(radar_baler.latitude['data']),  float(radar_baler.longitude['data'])),
    'SDAET': (float(radar_daet.latitude['data']),   float(radar_daet.longitude['data'])),
}
for name, (lat, lon) in radar_coords.items():
    ax.plot(lon, lat, marker='o', markersize=6, color='black',
            markerfacecolor='white', markeredgewidth=1.2, transform=proj, zorder=10)
    ax.text(lon + 0.15, lat + 0.15, name, fontsize=9, color='black',
            fontweight='bold', transform=proj,
            bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor='none'))

# 标题
ax.set_title('Radar Reflectivity Mosaic (2 km CAPPI)', fontsize=13, fontweight='bold', pad=10)
ax.text(0.5, -0.12, 'Longitude (E)', transform=ax.transAxes, ha='center', fontsize=11)
ax.text(-0.1, 0.5, 'Latitude (N)', transform=ax.transAxes, va='center', rotation='vertical', fontsize=11)

# 色条，放在图右侧，使用 inset_axes
cax = inset_axes(ax, width='3%', height='100%', loc='lower left',
                 bbox_to_anchor=(1.02, 0, 1, 1), bbox_transform=ax.transAxes)
cbar = plt.colorbar(mesh, cax=cax, ticks=level)
cbar.set_label('Reflectivity (dBZ)', fontsize=11, rotation=270, labelpad=18)
cbar.ax.tick_params(labelsize=9)

# 精确调整边距
fig.subplots_adjust(left=0.1, right=0.88, bottom=0.08, top=0.92)
plt.savefig('radar_mosaic_pyart_cartopy.png', dpi=300, bbox_inches='tight', pad_inches=0.1)
plt.show()
print('图片已保存， radar_mosaic_pyart_cartopy.png (300 dpi)')
image
image

image

3. 完整代码

上面的演示用了 300×300 的网格，分辨率够看效果但还是糙了点。下面这段是把 grid_shape 提到 800×800、并把绘图扩到吕宋岛全岛范围的版本，可直接复制运行。

代码语言：javascript
复制
# =====================================================
# Py-ART 雷达组网完整代码（Cartopy 出版级版本）
# =====================================================
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import cartopy.crs as ccrs
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import pyart
import cmaps

# 1. 读取雷达数据
radar_baguio = pyart.io.read_cfradial('20251103051105_SBAG_CAPPI_2km.nc')
radar_baler  = pyart.io.read_cfradial('20251103051001_SBAL_CAPPI_2km.nc')
radar_daet   = pyart.io.read_cfradial('20251103051038_SDAET_CAPPI_2km.nc')

# 2. 创建 GateFilter（可选）
def make_filter(radar):
    gf = pyart.filters.GateFilter(radar)
    gf.exclude_transition()
    gf.exclude_above('Reflectivity_CAPPI_2km', 100)
    return gf

# 3. 网格化组网
grid = pyart.map.grid_from_radars(
    (radar_daet, radar_baguio, radar_baler),
    grid_shape=(1, 800, 800),
    grid_limits=((0, 2000), (-2000000, 2000000), (-1000000, 1000000)),
    fields=['Reflectivity_CAPPI_2km'],
    weighting_function='barnes2',
    gridding_algo='map_gates_to_grid',
    roi_func='dist_beam',
)

# 4. 提取数据
refl = grid.fields['Reflectivity_CAPPI_2km']['data'][0]
lats = grid.point_latitude['data'][0]
lons = grid.point_longitude['data'][0]

# 5. 使用 Cartopy 绘图（出版级）
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['Times New Roman', 'DejaVu Serif']
plt.rcParams['mathtext.fontset'] = 'stix'
plt.rcParams['font.size'] = 11

cmap = cmaps.radar_1
bounds = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
norm = mcolors.BoundaryNorm(bounds, cmap.N, extend='both')

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 8), dpi=300)
ax = plt.axes(projection=proj)
ax.set_extent([116, 128, 12, 20], crs=proj)

mesh = ax.pcolormesh(lons, lats, refl, transform=proj, cmap=cmap, norm=norm, shading='auto')
ax.coastlines(resolution='10m', linewidth=0.6, color='dimgray')

gl = ax.gridlines(
    draw_labels=True, linewidth=0.3, color='gray', alpha=0.4, linestyle='--',
    xlocs=range(116, 129, 2), ylocs=range(12, 21, 2),
)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}

for name, radar in [('SBAG', radar_baguio), ('SBAL', radar_baler), ('SDAET', radar_daet)]:
    lat = float(radar.latitude['data'])
    lon = float(radar.longitude['data'])
    ax.plot(lon, lat, marker='o', markersize=6, color='black',
            markerfacecolor='white', markeredgewidth=1.2, transform=proj, zorder=10)
    ax.text(lon + 0.15, lat + 0.15, name, fontsize=9, fontweight='bold', transform=proj,
            bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor='none'))

ax.set_title('Radar Reflectivity Mosaic (2 km CAPPI)', fontsize=13, fontweight='bold', pad=10)
ax.text(0.5, -0.12, 'Longitude (E)', transform=ax.transAxes, ha='center', fontsize=11)
ax.text(-0.1, 0.5, 'Latitude (N)', transform=ax.transAxes, va='center', rotation='vertical', fontsize=11)

cax = inset_axes(ax, width='2.5%', height='100%', loc='lower left',
                 bbox_to_anchor=(1.02, 0, 1, 1), bbox_transform=ax.transAxes)
cbar = plt.colorbar(mesh, cax=cax, ticks=bounds)
cbar.set_label('Reflectivity (dBZ)', fontsize=11, rotation=270, labelpad=18)
cbar.ax.tick_params(labelsize=9)

fig.subplots_adjust(left=0.1, right=0.88, bottom=0.08, top=0.92)
plt.savefig('radar_mosaic_pyart_cartopy.png', dpi=300, bbox_inches='tight', pad_inches=0.1)
plt.show()

附录，常见问题

Q1，组网结果出现大片空白区域

先检查 grid_limits 是否覆盖了所有雷达的探测范围。雷达阵列中心是参考原点，三部雷达的相对位置如果偏到一边，默认 ±1000 km 的范围可能正好把某部雷达切到边缘。

另外 grid_shape 不宜过小，500 km 范围内 100×100 的网格分辨率太低，插值结果会出现明显的网格化条纹。

Q2，如何提高计算速度

三个方向。一是把 grid_shape 从 800×800 降到 400×400，本身就是 4 倍速度的差距；二是 roi_func 用 'constant' 代替 'dist_beam'，省掉每个网格点的距离计算；三是上 Dask 并行。

代码语言：javascript
复制
from dask.distributed import Client, LocalCluster
cluster = LocalCluster()
client = Client(cluster)

Q3，Py-ART 适合中国雷达数据吗

中国新一代天气雷达（CINRAD）数据通常是 SA/SB/CB 格式的二进制文件，Py-ART 原生不支持。可以先用 pycwr 或 pycinrad 读取 + 转格式（如 CF/Radial 或 ODIM_H5），再交给 Py-ART 处理。这条链路目前是国内最常见的方案。

参考资料

  • Py-ART 官方文档， https://arm-doe.github.io/pyart/
  • grid_from_radars API， https://arm-doe.github.io/pyart/API/generated/pyart.map.grid_from_radars.html
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2026-05-09，如有侵权请联系 cloudcommunity@tencent.com 删除
grid
对象
函数
教程
数据

本文分享自 气python风雨 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ，欢迎热爱写作的你一起参与！

grid
对象
函数
教程
数据
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
目录
  • 雷达系列 | Py-ART 多雷达组网实战（grid_from_radars）
    • 阅读路径
    • 1. 数据介绍
    • 2. Py-ART 组网
      • 2.1 读取多部雷达数据
      • 2.2 数据质控（GateFilter）
      • 2.3 网格化与组网（grid_from_radars）
      • 2.4 提取并绘制组网结果（Cartopy）
    • 3. 完整代码
    • 附录，常见问题
      • Q1，组网结果出现大片空白区域
      • Q2，如何提高计算速度
      • Q3，Py-ART 适合中国雷达数据吗
    • 参考资料
领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 粤公网安备44030502008569号

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

Copyright © 2013 - 2026 Tencent Cloud.

All Rights Reserved. 腾讯云 版权所有

登录 后参与评论