Matplotlib中混合Sigma坐标下非统一Y轴等高线绘制及地形遮挡处理

问题背景 #

需要在Matplotlib中绘制混合Sigma压力坐标的纬度-压力网格经向风场：

• 每个纬度（lat_plot，36×1数组）对应独立的压力层数组（level_sigma，36×26数组）
• 要求等高线从地形上方开始绘制，避免与地形区域重叠
• 现有代码仅支持统一Y轴压力层的等高线绘制，无法适配混合Sigma坐标需求

实现方案 #

1. 地形压力转换与数据掩码 #

首先将地形高度转换为对应压力值（示例中简化处理，实际需根据气象模型静力学公式调整），然后对经向风数据进行掩码，把地形下方的数值设为np.nan，Matplotlib会自动跳过这些区域的绘制。

2. 构建可变Y轴的网格坐标 #

混合Sigma坐标下每个纬度的压力层不同，需构建对应维度的X、Y网格：

• X网格：基于lat_plot生成与level_sigma同维度的纬度网格
• Y网格：直接使用level_sigma作为每个纬度对应的压力层网格

3. 适配可变Y轴的绘图函数 #

向contourf/contour传入二维的X、Y网格（而非一维数组），让Matplotlib识别每个纬度对应的独立压力层。

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 生成基础数据
lat_plot = np.linspace(-90, -60, 36)
level = np.linspace(0, 1000, 26)[::-1]
v = np.random.rand(len(lat_plot), len(level))
topo_zm = np.random.uniform(0, 3000, len(lat_plot))

# 生成混合Sigma坐标的压力层（模拟数据）
level_sigma = np.zeros_like(v)
for i in range(len(level)):
    a = np.random.randint(0, 1000, len(level_sigma))[::-1]
    order = a.argsort()
    level_sigma[:,i] = a[order]

# 地形压力转换（示例简化：将地形高度转为hPa量级，实际需用静力学公式计算）
topo_pressure = topo_zm / 10

# 数据掩码：将地形下方的风场数据设为NaN，避免绘制地形区域的等高线
for i in range(len(lat_plot)):
    below_topo = level_sigma[i, :] < topo_pressure[i]
    v[i, below_topo] = np.nan

# 绘图配置
levels = [-1, -0.875, -0.75, -0.625, -0.5, -0.375, -0.25, -0.125, 0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1]
ticks = [-1, -0.5, 0, 0.5, 1]

fig = plt.figure(figsize=[15,7.5])
ax = fig.add_subplot()

# 构建匹配维度的X、Y网格（indexing='ij'确保与lat_plot×level维度一致）
X, _ = np.meshgrid(lat_plot, level, indexing='ij')
Y = level_sigma

# 绘制填色等高线和轮廓线
im2 = ax.contourf(X, Y, v, levels=levels, cmap="coolwarm", extend="both")
im1 = ax.contour(X, Y, v, levels=levels, colors="black", linewidths=0.5, linestyles="solid")

# 坐标轴配置
ax.set_xlim(-90, -60)
ax.set_xlabel("latitude [°]")
ax.set_xticks([-90, -80, -70, -60])
ax.invert_yaxis()
ax.set_ylim(1000, 300)
ax.set_ylabel("pressure [hPa]")
ax.set_yticks([900, 700, 500, 300])
ax.set_title("4PIC", fontsize=10)

# 绘制地形阴影
ax1 = ax.twinx()
ax1.plot(lat_plot, topo_zm, color="Black")
ax1.fill_between(lat_plot, topo_zm, color="Black", alpha=0.3)
ax1.set_ylim(0,10000)
ax1.axes.get_yaxis().set_visible(False)

fig.suptitle("Zonal-mean meriodional summer (DJ) wind averaged around West- (180°-360°E) and East-Antarctica (0°-180°E)", x=0.52, y=0.94)

# 添加色标
cbar = fig.colorbar(im2, ax=ax, ticks=ticks)
cbar.set_label("meridional wind [m/s]")

plt.show()

关键说明 #

• 数据掩码是核心：通过np.nan标记地形下方的无效数据，确保等高线仅在地形上方绘制
• 二维网格传递：向绘图函数传入匹配维度的X、Y网格，适配混合Sigma坐标的可变压力层
• 地形压力转换：示例为简化处理，实际应用需根据气象模型的静力学关系计算准确的地形对应压力值

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