NetCDF数据4D曲面绘图问题：仅显示单个层级如何修正？

我有一个包含time、lat、lon、lev四个维度、5个变量（DU001至DU005）的NetCDF数据，想要绘制4D曲面图：以lat为x轴、lon为y轴、lev为z轴，用颜色表示变量DU001（共72个层级）。但运行代码后只显示单个曲面，推测仅读取了一个层级，求修正方法。

错误代码如下：

import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import cm
path= 'D:\\DATA\\2015\\test_data'# Open the NetCDF file
data = xr.open_dataset('D:\\DATA\\2015\\test_data\\MERRA2_400.inst3_3d_aer_Nv.20150515.SUB.nc')
# Select the DU01 variable and the lat, long, and lev dimensions
lat = data['lat']
lon = data['lon']
lev = data['lev']
DMR = data['DU001'] 
# Reshape the data
du001_2d = DMR[:, :, :].squeeze()
dmr_values = du001_2d.values.squeeze()
# Create meshgrid for coordinates
lon_2d, lat_2d = np.meshgrid(lon, lat)
# Create the 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the surface for each level
for i, level in enumerate(lev):
    ax.plot_surface(lon_2d, lat_2d, dmr_values[i], cmap='viridis')
    # Set labels and title
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    ax.set_zlabel('Level')
    ax.set_title('3D Plot of Dust Mixing Ratio')

    # Set the z-limits based on the valid range of the lev array
    ax.set_zlim(lev[0], lev[71])  # Assuming lev is a 1D array

# Display the plot
plt.show()

问题根源 #

1. z轴赋值错误：你把DU001的数值当作z轴数据传入plot_surface，但z轴应该对应lev的层级高度，变量值要通过颜色来体现。
2. 维度处理不当：无差别使用squeeze()可能误删必要维度；循环内重复设置标签、标题和z轴范围，完全冗余。
3. 曲面遮挡问题：默认曲面不透明，后面的层级会被前面的挡住，无法看到多层效果。

修正后的代码 #

import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import cm

# 读取NetCDF数据
data = xr.open_dataset('D:\\DATA\\2015\\test_data\\MERRA2_400.inst3_3d_aer_Nv.20150515.SUB.nc')

# 提取变量和维度
lat = data['lat'].values
lon = data['lon'].values
lev = data['lev'].values
# 明确去掉time维度（假设仅单个时间点）
dmr_values = data['DU001'].squeeze(dim='time').values

# 创建经纬度网格
lon_2d, lat_2d = np.meshgrid(lon, lat)

# 创建3D绘图对象
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 归一化函数：将数据映射到0-1范围，适配颜色板
def normalize(data):
    return (data - data.min()) / (data.max() - data.min())

# 遍历每个层级绘制曲面
for i in range(len(lev)):
    # z轴固定为当前层级高度，变量值通过facecolors传递为颜色
    surf = ax.plot_surface(lon_2d, lat_2d, np.full_like(lat_2d, lev[i]), 
                           facecolors=cm.viridis(normalize(dmr_values[i])),
                           alpha=0.5,  # 设置半透明，避免遮挡下层
                           rstride=1, cstride=1)

# 设置坐标轴标签和标题
ax.set_xlabel('经度')
ax.set_ylabel('纬度')
ax.set_zlabel('层级')
ax.set_title('DU001 多层3D分布')

# 设置z轴范围匹配层级区间
ax.set_zlim(lev.min(), lev.max())

# 添加颜色条，关联DU001的数值范围
norm = plt.Normalize(dmr_values.min(), dmr_values.max())
sm = cm.ScalarMappable(norm=norm, cmap='viridis')
sm.set_array([])
fig.colorbar(sm, ax=ax, label='DU001 数值')

plt.show()

关键修正说明 #

• z轴数据修正：用np.full_like(lat_2d, lev[i])生成和经纬度网格同形状的数组，每个层级的z轴固定为当前lev值，变量值通过facecolors映射为颜色。
• 透明度设置：添加alpha=0.5让曲面半透明，能同时展示多层分布。
• 颜色映射优化：通过归一化函数将DU001数值适配到颜色板的0-1区间，添加颜色条明确数值与颜色的对应关系。
• 维度清理：明确指定squeeze(dim='time')去掉时间维度，避免误删其他必要维度。
• 代码效率优化：把标签、标题、z轴范围的设置移到循环外，只执行一次。

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