雷达系列 | wradlib 多雷达组网实战

!pip install wradlib xradar --upgrade -i https://pypi.mirrors.ustc.edu.cn/simple/

# 导入常用库
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import wradlib as wrl
import xradar as xd
import pyproj
import cmaps
import warnings
warnings.filterwarnings('ignore')

import cartopy.crs as ccrs
import cartopy.feature as cfeature

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

print('依赖检测完成')

# 步骤 1，构造 UTM 投影器
proj_utm = pyproj.CRS.from_epsg(32651)                                          # UTM Zone 51N
transformer = pyproj.Transformer.from_crs('EPSG:4326', proj_utm, always_xy=True)

# 步骤 2，逐个雷达读取 + 构造带全局 UTM 坐标的 DataArray
radar_datasets = {}
for filepath, name in [(file_baguio, 'SBAG'), (file_baler, 'SBAL'), (file_daet, 'SDAET')]:
    dtree = xd.io.open_cfradial1_datatree(filepath)
    dtree_geo = xd.georeference.get_x_y_z_tree(dtree)
    sweep = dtree_geo['/sweep_0']

    # 本地切平面坐标，原点在雷达站
    x_local = sweep['x'].values
    y_local = sweep['y'].values
    refl = sweep['Reflectivity_CAPPI_2km']

    # 雷达站经纬度 -> UTM
    radar_lon = float(sweep.coords['longitude'].values)
    radar_lat = float(sweep.coords['latitude'].values)
    radar_x, radar_y = transformer.transform(radar_lon, radar_lat)

    # 创建带全局 UTM 坐标的 DataArray
    da = xr.DataArray(
        refl.values,
        dims=['azimuth', 'range'],
        coords={
            'x': (['azimuth', 'range'], x_local + radar_x),
            'y': (['azimuth', 'range'], y_local + radar_y),
        },
        name='Reflectivity',
    )
    # 过滤 fill_value（原始数据里超大值都是无效）
    da = da.where(da < 1000)

    radar_datasets[name] = {
        'da': da,
        'site': (radar_x, radar_y),
        'site_lonlat': (radar_lon, radar_lat),
    }
    print(f'{name}, x={float(da.x.min().values):.0f}~{float(da.x.max().values):.0f}, '
          f'y={float(da.y.min().values):.0f}~{float(da.y.max().values):.0f}')

# 步骤 3，从所有雷达数据的 bbox 推导公共网格范围
xmin = min([rd['da'].x.min().values for rd in radar_datasets.values()])
xmax = max([rd['da'].x.max().values for rd in radar_datasets.values()])
ymin = min([rd['da'].y.min().values for rd in radar_datasets.values()])
ymax = max([rd['da'].y.max().values for rd in radar_datasets.values()])
print(f'数据 bbox， x={xmin:.0f}~{xmax:.0f}, y={ymin:.0f}~{ymax:.0f}')

# 目标网格 300x300，覆盖全部数据
x = np.linspace(xmin, xmax, 300)
y = np.linspace(ymin, ymax, 300)
cart = xr.Dataset(coords={'x': (['x'], x), 'y': (['y'], y)})
print(f'公共网格， {len(x)} x {len(y)}')

# 顺手算一下经纬度网格，绘图时直接用
xx_grid, yy_grid = np.meshgrid(x, y)
transformer_inv = pyproj.Transformer.from_crs(proj_utm, 'EPSG:4326', always_xy=True)
lons_grid, lats_grid = transformer_inv.transform(xx_grid, yy_grid)

# pulse_volume 参数
#   range_resolution    距离分辨率（米），CF/Radial 元数据里读出来一般是 250~1000 m
#   beamwidth           波束宽度（度），常见 S 波段雷达约 1.0°
RANGE_RES = 504.0
BEAM_WIDTH = 1.0

radar_grids = []
quality_grids = []

for name, rd in radar_datasets.items():
    da = rd['da']
    site = rd['site']

    # 采样体积作为质量因子（体积越小越可靠）
    pv = da.wrl.qual.pulse_volume(RANGE_RES, BEAM_WIDTH)

    # 探测半径 = 距雷达站最远像素的距离
    max_range = float(np.nanmax(np.sqrt((da.x - site[0]) ** 2 + (da.y - site[1]) ** 2)))
    # 用数据中心做插值参考中心，注意 togrid 期望 (y, x) 顺序
    center = (float(np.nanmean(da.y.values)), float(np.nanmean(da.x.values)))

    # 插值到公共网格（IDW）
    g = da.wrl.comp.togrid(cart, radius=max_range, center=center, interpol=wrl.ipol.Idw)
    q = pv.wrl.comp.togrid(cart, radius=max_range, center=center, interpol=wrl.ipol.Idw)

    radar_grids.append(g)
    quality_grids.append(1.0 / (q + 0.001))                         # 体积越小，质量越高

    print(f'{name}, gridded={float(np.nanmin(g.values)):.2f}~{float(np.nanmax(g.values)):.2f} dBZ')

print('所有雷达网格化完成')

# 沿新维度 'radar' 堆叠
radar = xr.DataArray(list(radar_datasets.keys()), dims='radar')
radargrids = xr.concat(radar_grids, dim=radar)
qualitygrids = xr.concat(quality_grids, dim=radar)

# Knockout 组合
composite_ko = radargrids.wrl.comp.compose_ko(qualitygrids)
# 加权组合
composite_weighted = radargrids.wrl.comp.compose_weighted(qualitygrids)

print(f'Knockout 组合范围， {float(np.nanmin(composite_ko.values)):.2f} ~ {float(np.nanmax(composite_ko.values)):.2f} dBZ')
print(f'加权组合范围，   {float(np.nanmin(composite_weighted.values)):.2f} ~ {float(np.nanmax(composite_weighted.values)):.2f} dBZ')

# 步骤 5，将 UTM 网格转换回经纬度并绘图
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

fig, axes = plt.subplots(
    1, 2, figsize=(16, 7), dpi=300,
    subplot_kw={'projection': ccrs.PlateCarree()},
)

for idx, (ax, comp_data_raw, title) in enumerate(zip(
    axes,
    [composite_ko, composite_weighted],
    ['(a) Knockout composition', '(b) Weighted composition'],
)):
    ax.set_extent([120, 126, 12, 16.5], crs=ccrs.PlateCarree())

    # 把无效值（<=0）置 nan，避免显示为黑底
    comp_data = np.where(comp_data_raw.values > 0, comp_data_raw.values, np.nan)

    mesh = ax.contourf(
        lons_grid, lats_grid, comp_data,
        transform=ccrs.PlateCarree(), 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'}
    if idx == 1:                        # 右图不重复 y 轴标签
        gl.left_labels = False

    # 雷达站位置直接用 radar_datasets 里记录的经纬度
    for r_name, rd in radar_datasets.items():
        lon, lat = rd['site_lonlat']
        ax.plot(lon, lat, marker='o', markersize=5, color='black',
                markerfacecolor='white', markeredgewidth=1.0,
                transform=ccrs.PlateCarree(), zorder=10)
        ax.text(lon + 0.15, lat + 0.15, r_name, fontsize=8, color='black',
                fontweight='bold', transform=ccrs.PlateCarree(),
                bbox=dict(boxstyle='round,pad=0.15', facecolor='white',
                          alpha=0.7, edgecolor='none'))

    ax.set_title(title, fontsize=12, fontweight='bold', pad=8)

# 共享色条，放在两个子图下方中央
cbar_ax = fig.add_axes([0.35, 0.06, 0.3, 0.025])
cbar = fig.colorbar(mesh, cax=cbar_ax, orientation='horizontal', ticks=level)
cbar.set_label('Reflectivity (dBZ)', fontsize=11)
cbar.ax.tick_params(labelsize=9)

fig.suptitle('wradlib Radar Reflectivity Mosaic (2 km CAPPI)',
             fontsize=13, fontweight='bold', y=0.98)

fig.subplots_adjust(left=0.06, right=0.94, bottom=0.14, top=0.90, wspace=0.12)
plt.savefig('radar_mosaic_wradlib_cartopy.png', dpi=300, bbox_inches='tight', pad_inches=0.1)
plt.show()
print('图片已保存， radar_mosaic_wradlib_cartopy.png (300 dpi)')