如何用无经纬度5km格点数据结合Cartopy绘制雷达反射率子图？

核心问题分析 #

你的代码无法显示反射率数据，主要有三个关键问题：

• 子图未使用Cartopy投影轴：默认创建的matplotlib轴无法正确解析地理坐标转换参数，导致反射率数据无法映射到正确位置。
• 经纬度维度不匹配：pcolormesh需要2D网格数据，但你转换后的经纬度是1D数组，无法和2D的反射率数据匹配。
• 经纬度数据来源错误：你使用了外部的xg['x']/xg['y']转换经纬度，而非当前读取文件的x/y值，可能导致坐标和数据不对应。

修正步骤 #

1. 创建Cartopy投影子图：在plt.subplots中通过subplot_kw指定投影类型，确保每个子图都是GeoAxes。
2. 生成2D经纬度网格：对每个文件的x、y值做meshgrid转换，得到和反射率维度一致的2D经纬度数组。
3. 在循环内处理坐标转换：针对每个读取的文件，单独提取x、y并转换为经纬度，保证数据对应。
4. 设置轴范围：根据转换后的经纬度范围限制显示区域，确保数据在可视范围内。

修正后的完整代码 #

import pyproj
import numpy as np
import pyart
import matplotlib.pyplot as plt
import glob
import cartopy.crs as ccrs
import cartopy.feature as feat

center_lat = 26.820560455322266
center_lon = 75.81861114501953
# 定义AEQD投影（雷达站中心的方位等距投影）
proj = pyproj.Proj(proj='aeqd', lat_0=center_lat, lon_0=center_lon)

cff = glob.glob("C:/Users/shrey/RADAR_DATA/RADAR_DATA_JAIPUR/JPR220524IMD-B/output/outgrid_new/grid_new*nc")
plot_proj = ccrs.PlateCarree()  # 绘图使用的投影

for files in cff[:]:
    read_files = pyart.io.read_grid(files)
    xg_files = read_files.to_xarray()
    
    # 从当前文件获取x、y值（1D数组）
    x_values = xg_files['x'].values
    y_values = xg_files['y'].values
    
    # 生成2D网格的x、y数组
    x_grid, y_grid = np.meshgrid(x_values, y_values)
    # 转换为2D的经纬度数组
    lon_grid, lat_grid = proj(x_grid, y_grid, inverse=True)
    
    # 计算显示范围（Cartopy的extent格式为[lon_min, lon_max, lat_min, lat_max]）
    lat_min, lat_max = np.min(lat_grid), np.max(lat_grid)
    lon_min, lon_max = np.min(lon_grid), np.max(lon_grid)
    lat_lon_extent = [lon_min, lon_max, lat_min, lat_max]
    
    z_heights = [500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000]
    num_subplots = len(z_heights)
    num_rows = num_subplots // 2
    num_cols = 2

    # 创建带Cartopy投影的子图
    fig, axes = plt.subplots(
        nrows=num_rows, 
        ncols=num_cols, 
        figsize=(12, 12),
        subplot_kw={'projection': plot_proj}
    )
    axes = axes.flatten()

    for i, z_height in enumerate(z_heights):
        ax = axes[i]
        # 选择对应高度的反射率数据
        reflectivity = xg_files['REF'].sel(z=z_height, method='nearest').values
        
        # 绘制反射率数据，指定数据的坐标投影
        im = ax.pcolormesh(
            lon_grid, 
            lat_grid, 
            reflectivity, 
            vmin=0, 
            vmax=60, 
            cmap='rainbow', 
            transform=plot_proj
        )
        
        # 设置轴显示范围
        ax.set_extent(lat_lon_extent, crs=plot_proj)
        
        # 添加地理要素
        ax.add_feature(feat.STATES, linewidth=0.5)
        ax.add_feature(feat.BORDERS, linewidth=1)
        ax.set_title(f"Reflectivity at z = {z_height} m")

    fig.suptitle(f"File: {files.split('_')[-1]}", fontsize=16)
    plt.tight_layout()

    # 添加色标
    cbar = fig.colorbar(im, ax=axes, orientation='vertical', shrink=0.8)
    cbar.set_label('Reflectivity (dBZ)')

    plt.show()

额外说明 #

• Cartopy的set_extent参数顺序为[lon_min, lon_max, lat_min, lat_max]，需注意和你之前的格式区分。
• 每个文件单独处理坐标转换，确保x/y和反射率数据来自同一文件，避免数据不匹配。
• 使用meshgrid生成2D网格是核心修正点，这样经纬度数组的维度才能和反射率的(y,x)维度完全对应。

内容的提问来源于stack exchange，提问作者Shreyasi Upadhyay