雷达系列 | 两种雷达坐标转换经纬度方法对比分析

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# MetPy相关
from metpy.calc import azimuth_range_to_lat_lon
from metpy.io import Level3File
from metpy.units import units

# wradlib相关
from wradlib.georef import polar

print("所有库导入成功!")

所有库导入成功!

# 读取真实MetPy雷达数据
filename = 'KOUN_SDUS54_N0QTLX_201305202016 (1)'
f = Level3File(filename)

# 从文件中提取数据
datadict = f.sym_block[0][0]
data = f.map_data(datadict['data'])

# 获取方位角和距离
az = units.Quantity(np.array(datadict['start_az'] + [datadict['end_az'][-1]]), 'degrees')
rng = units.Quantity(np.linspace(0, f.max_range, data.shape[-1] + 1), 'kilometers')

# 获取雷达站点位置
center_lon = f.lon
center_lat = f.lat

print(f"数据文件: {filename}")
print(f"数据形状: {data.shape}")
print(f"方位角数量: {len(az)}")
print(f"距离库数量: {len(rng)} (范围: 0-{f.max_range} km)")
print(f"雷达中心位置: 经度={center_lon:.4f}°, 纬度={center_lat:.4f}°")

# 使用MetPy进行坐标转换
xlocs, ylocs = azimuth_range_to_lat_lon(az, rng, center_lon, center_lat)
print('经度范围：{:.3f} °E  ~  {:.3f} °E'.format(xlocs.min(), xlocs.max()))
print('纬度范围：{:.3f} °N  ~  {:.3f} °N'.format(ylocs.min(), ylocs.max()))

数据文件: KOUN_SDUS54_N0QTLX_201305202016 (1)
数据形状: (360, 460)
方位角数量: 361
距离库数量: 461 (范围: 0-460.0 km)
雷达中心位置: 经度=-97.2780°, 纬度=35.3330°
经度范围：-102.351 °E  ~  -92.205 °E
纬度范围：31.196 °N  ~  39.470 °N

# wradlib方法：三维球面坐标转换
# 准备与MetPy匹配的坐标数据
cent_lon = f.lon
cent_lat = f.lat
# 雷达站点坐标 (经度, 纬度, 高度)
csite=(cent_lon,cent_lat,0)
theta= np.full(1, 361)
# 创建网格坐标
coords = polar.spherical_to_proj(rng*1000, az, theta, csite)
print(f"wradlib转换完成")
print(f"转换后坐标形状: {coords.shape}")
print(f"X坐标范围: {coords[..., 0].min():.2f} - {coords[..., 0].max():.2f}")
print(f"Y坐标范围: {coords[..., 1].min():.2f} - {coords[..., 1].max():.2f}")
print(f"Z坐标范围: {coords[..., 2].min():.2f} - {coords[..., 2].max():.2f}")

wradlib转换完成
转换后坐标形状: (361, 461, 3)
X坐标范围: -102.33 - -92.23
Y坐标范围: 31.19 - 39.47
Z坐标范围: 0.00 - 20458.30

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.plots import colortables, add_timestamp
import cartopy
# ---------- 色表与投影 ----------
ctable = ('NWSStormClearReflectivity', -20, 0.5)   # dBZ
norm, cmap = colortables.get_with_steps(*ctable)
crs_map = ccrs.LambertConformal()                  # 地图投影

# ---------- 画布 ----------
fig = plt.figure(figsize=(15, 8))                  # 左右 2 子图，横向排布

# ---------- 通用装饰 ----------
def draw_one(ax, x, y, data, title):
    """画一张图并统一装饰"""
    mesh = ax.pcolormesh(x, y, data, norm=norm, cmap=cmap,shading='auto', 
                         transform=ccrs.PlateCarree())
    ax.set_extent([cent_lon-3, cent_lon+3,
                   cent_lat-3, cent_lat+3])
    ax.set_aspect('equal', 'datalim')
    ax.set_title(title, fontsize=14)

    # 经纬度网格 / 刻度
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                      linewidth=1, color='gray', alpha=0.5, linestyle='--',
                      x_inline=False, y_inline=False, rotate_labels=False)
    gl.top_labels = gl.right_labels = False
    gl.xformatter = cartopy.mpl.gridliner.LONGITUDE_FORMATTER
    gl.yformatter = cartopy.mpl.gridliner.LATITUDE_FORMATTER
    gl.xlabel_style = gl.ylabel_style = {'size': 11}

    # 雷达站位置
    ax.plot(cent_lon, cent_lat, 'r+', markersize=12, label='雷达站')
    ax.legend(loc='upper right')

    # 时间戳（用 MetPy 自带函数）
    add_timestamp(ax, f.metadata['prod_time'], y=0.02, high_contrast=True)

    return mesh

# ---------- 左侧：MetPy ----------
ax1 = fig.add_subplot(1, 2, 1, projection=crs_map)
m1 = draw_one(ax1, xlocs, ylocs, data, 'MetPy 方法')

# ---------- 右侧：wradlib ----------
ax2 = fig.add_subplot(1, 2, 2, projection=crs_map)
# 如果 wradlib 的坐标维度是 (R, A, 2)
plot_coords = coords 
lons = plot_coords[:, :, 0]
lats = plot_coords[:, :, 1]

m2 = draw_one(ax2, lons, lats, data, 'wradlib 方法')

# ---------- 统一 colorbar ----------
cbar = plt.colorbar(m1, ax=[ax1, ax2], orientation='horizontal',
                    pad=0.06, fraction=0.05)
cbar.set_label('Reflectivity (dBZ)', fontsize=12)
plt.tight_layout()
plt.show()

/tmp/ipykernel_167/1266596128.py:61: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()