雷达系列 | 如何转换雷达三维球面坐标为经纬度坐标系

import numpy as np
from wradlib.georef import polar

# 定义距离、方位角和仰角
r = np.array([0., 0., 111., 111., 111., 111.]) * 1000
az = np.array([0., 180., 0., 90., 180., 270.])
th = np.array([0., 0., 0., 0., 0., 0.5])

# 定义站点位置
csite = (9.0, 48.0, 0)

# 执行坐标转换
coords = polar.spherical_to_proj(r, az, th, csite)

# 输出结果
for coord in coords:
    print('{0:7.4f}, {1:7.4f}, {2:7.4f}'.format(*coord))

 9.0000, 48.0000,  0.0000
 9.0000, 48.0000,  0.0000
 9.0000, 48.9981, 725.7160
10.4872, 47.9904, 725.7160
 9.0000, 47.0017, 725.7160
 7.5131, 47.9904, 1694.2234

import cartopy.crs as ccrs
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np

from metpy.calc import azimuth_range_to_lat_lon
from metpy.cbook import get_test_data
from metpy.io import Level3File
from metpy.plots import add_timestamp, colortables, USCOUNTIES
from metpy.units import units

# 打开雷达数据文件
name = '/home/mw/project/KOUN_SDUS54_N0QTLX_201305202016'
f = Level3File(name)

# 从文件对象中提取数据字典
datadict = f.sym_block[0][0]

# 根据文件指定的比例尺将数据转换为数组
data = f.map_data(datadict['data'])
print("雷达原始数据：",data)
# 获取方位角和距离，并为其分配单位
az = units.Quantity(np.array(datadict['start_az'] + [datadict['end_az'][-1]]), 'degrees')
rng =   np.linspace(0, f.max_range, data.shape[-1] + 1)*1000
print("方位角",az)
print("距离",rng)
print(az.shape)
# 从文件中提取中心经度和纬度
cent_lon = f.lon
cent_lat = f.lat

##设置参数
csite=(cent_lon,cent_lat,0)
theta= np.full(3, 361)

coords = polar.spherical_to_proj(rng, az, th, csite)
coords.shape

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

# 假设我们需要绘制第一个维度的结果（根据您的输出形状 (6, 361, 461, 3)）
# 取出经度 (lon) 和 纬度 (lat)
# coords[..., 0] 是经度, coords[..., 1] 是纬度, coords[..., 2] 是高度
# 这里我们需要将数据降维以匹配 pcolormesh 的输入要求
# 假设 data 是 (361, 461) 的 2D 数组，我们取 coords 的第一个 slice
plot_coords = coords[0] 
lons = plot_coords[:, :, 0]
lats = plot_coords[:, :, 1]

# 确保 data 的维度与 lons/lats 匹配
# 注意：pcolormesh 的坐标网格通常比数据多一行/一列，或者维度相同（取决于 shading='nearest' 或 'flat'）
# 这里我们假设维度直接对应
plot_data = data  # 需确保 plot_data 是 2D 数组 (361, 461)

# 创建画布和投影
fig = plt.figure(figsize=(12, 10))
# 使用 PlateCarree 投影（等距圆柱投影），因为我们的 coords 是经纬度
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# 添加地图特征
ax.add_feature(cfeature.COASTLINE, linewidth=1)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.STATES, linestyle='--', linewidth=0.5) # 添加州界/省界
ax.gridlines(draw_labels=True, linestyle='--', alpha=0.5)

# 设置雷达回波的色标 (Reflectivity)
# 使用 MetPy 自带的色表，或者使用 standard 'pyart_HomeyerRainbow' 等
cmap = plt.get_cmap('gist_ncar') 
cmap.set_bad('white', alpha=0) # 将 NaN 值设为透明

# 绘制数据
# transform=ccrs.PlateCarree() 告诉 cartopy 输入的 lons 和 lats 是经纬度数据
mesh = ax.pcolormesh(lons, lats, plot_data, 
                     shading='auto', 
                     cmap=cmap, 
                     vmin=0, vmax=70, # 反射率因子通常在 0-70 dBZ 之间
                     transform=ccrs.PlateCarree())

# 添加色标
cbar = plt.colorbar(mesh, ax=ax, orientation='vertical', fraction=0.046, pad=0.04)
cbar.set_label('Reflectivity (dBZ)')

# 标记雷达站点位置
site_lon, site_lat, _ = csite
ax.plot(site_lon, site_lat, 'r+', markersize=15, transform=ccrs.PlateCarree(), label='Radar Site')
ax.legend()

# 设置标题
ax.set_title('KOUN Radar Reflectivity (Projected from Spherical Coords)', fontsize=14)

plt.show()