如何提取并绘制所有CMIP6模型均处于集合均值1标准差范围内的区域

看起来你已经搞定了CMIP6集合均值的地图绘制，现在想标记出所有模型都落在均值±1标准差范围内的区域——也就是模型输出变异最小、一致性最高的区域对吧？你的初步思路方向是对的，但用循环逐个处理模型效率太低，而且逻辑上有点小问题，咱们用numpy的向量化操作就能轻松搞定这个需求～

核心思路拆解 #

我们的目标是生成一个布尔掩码数组，标记所有模型都满足「值在集合均值±1标准差区间内」的网格点，然后把这个掩码区域用阴影叠加到现有地图上。

第一步：准备所有模型的二维数据 #

首先你需要把每个模型的(33,180,360)原始数据，对时间维度（axis=0）做平均，得到每个模型的(180,360)二维数组，再把9个模型的二维数组堆叠成一个(9,180,360)的大数组：

# 假设你已经整理好9个模型的文件路径列表
model_paths = [
    'path/to/model_1.nc',
    'path/to/model_2.nc',
    # ... 剩下7个模型的路径
]

all_models_2d = []
for path in model_paths:
    nc = Dataset(path, 'r')
    # 读取单个模型的时间序列数据
    model_data = nc.variables['tauu'][:]  # shape (33,180,360)
    # 对时间维度求平均，得到二维空间场
    model_2d = np.mean(model_data, axis=0)
    all_models_2d.append(model_2d)
# 转换为numpy数组，shape变为(9,180,360)
all_models_2d = np.array(all_models_2d)

第二步：计算集合均值与标准差 #

这部分你应该已经完成了，这里再明确下对应逻辑（如果已经有集合均值文件，直接用你原来的TS2_mean_hist即可）：

# 对模型维度求均值，得到集合均值场
ensemble_mean = np.mean(all_models_2d, axis=0)  # shape (180,360)
# 计算模型间的标准差
stan_dev = np.std(all_models_2d, axis=0)  # shape (180,360)

# 计算均值±1标准差的区间边界
sd_minus_1 = ensemble_mean - stan_dev
sd_plus_1 = ensemble_mean + stan_dev

第三步：生成一致性掩码 #

用numpy的向量化操作替代循环，一次性检查所有模型的每个网格点是否在区间内，再筛选出所有模型都满足条件的网格点：

# 检查每个模型的每个网格点是否在[sd_minus_1, sd_plus_1]内
in_range = (all_models_2d >= sd_minus_1) & (all_models_2d <= sd_plus_1)
# 只有当所有模型都满足条件时，该网格点才标记为True
consistent_mask = np.all(in_range, axis=0)  # shape (180,360)

这里consistent_mask中的True就是你需要标记的低变异区域。

第四步：在现有地图上叠加阴影 #

把掩码转换为Basemap的投影坐标，然后用pcolor或contourf绘制阴影区域，设置alpha=0让颜色透明，只保留阴影纹理：

# 把经纬度网格转换为Basemap投影坐标
x_mask, y_mask = m(lon, lat)
# 绘制阴影，hatch参数可选'///'、'xxx'、'+++'等纹理样式
m.pcolor(x_mask, y_mask, consistent_mask, hatch='///', alpha=0, edgecolor='gray', linewidth=0.1)

完整修改后的代码 #

把上述步骤整合到你原有的绘图代码中，最终代码如下：

from netCDF4 import Dataset
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import numpy as np
import os
os.environ["PROJ_LIB"] = "C:/Users/username/miniconda3/Library/share;" #fixr
from mpl_toolkits.basemap import Basemap
from pylab import *
import matplotlib as mpl
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap

# ---------------------新增：加载并处理所有模型数据---------------------
model_paths = [
    'path/to/model_1.nc',
    'path/to/model_2.nc',
    # 补充剩余7个模型的路径
]

all_models_2d = []
for path in model_paths:
    nc = Dataset(path, 'r')
    model_data = nc.variables['tauu'][:]  # shape (33,180,360)
    model_2d = np.mean(model_data, axis=0)  # 时间维度平均
    all_models_2d.append(model_2d)
all_models_2d = np.array(all_models_2d)  # shape (9,180,360)

# 计算集合均值与标准差（如果已有集合均值文件，可直接用TS2_mean_hist替代ensemble_mean）
ensemble_mean = np.mean(all_models_2d, axis=0)
stan_dev = np.std(all_models_2d, axis=0)
sd_minus_1 = ensemble_mean - stan_dev
sd_plus_1 = ensemble_mean + stan_dev

# 生成一致性掩码
in_range = (all_models_2d >= sd_minus_1) & (all_models_2d <= sd_plus_1)
consistent_mask = np.all(in_range, axis=0)
# ---------------------新增部分结束---------------------

# 原有绘图代码
fig_index=1
fig = plt.figure(num=fig_index, figsize=(12,7), facecolor='w')
fbot_levels = arange(-0.3, 0.5,0.05)
fname_hist_tauu='C:/Users/userbame/Historical data analysis/Historical/Wind/tauu_hist_ensmean_so.nc'
ncfile_tauu_hist = Dataset(fname_hist_tauu, 'r', format='NETCDF4')
TS2_hist=ncfile_tauu_hist.variables['tauu'][:]
TS2_mean_hist = np.mean(TS2_hist, axis=(0))
LON=ncfile_tauu_hist.variables['LON'][:]
LAT=ncfile_tauu_hist.variables['LAT'][:]
ncfile_tauu_hist.close()
lon,lat=np.meshgrid(LON,LAT)
ax1 = plt.axes([0.1, 0.225, 0.5, 0.6])
meridians=[0,1,1,1]
m = Basemap(projection='spstere',lon_0=0,boundinglat=-35)
m.drawcoastlines()
x, y =m(lon,lat)
m.contourf(x,y,TS2_mean_hist , fbot_levels, origin='lower', cmap=cm.RdYlGn)
m.drawparallels(np.arange(-90.,120.,10.),labels=[1,0,0,0]) # draw parallels
m.drawmeridians(np.arange(0.,420.,30.),labels=meridians) # draw meridians

# ---------------------新增：叠加一致性区域阴影---------------------
x_mask, y_mask = m(lon, lat)
# 调整hatch参数可更换阴影样式，比如'xxx'、'+++'
m.pcolor(x_mask, y_mask, consistent_mask, hatch='///', alpha=0, edgecolor='gray', linewidth=0.1)
# ---------------------新增结束---------------------

coloraxis = [0.1, 0.1, 0.5, 0.035]
cx = fig.add_axes(coloraxis, label='m', title='Wind Stress/ Pa')
cbar=plt.colorbar(cax=cx,orientation='horizontal',ticks=list(fbot_levels))
plt.savefig('C:/Users/username/Historical data analysis/Historical/Wind/Wind_hist_AP.png')

小提示 #

• 如果你已经有集合均值的NC文件，直接用TS2_mean_hist替代ensemble_mean即可，不用重复计算；
• 阴影纹理可通过hatch参数调整，比如'\\\\\\'（反斜线）、'....'（点）等，选你觉得最清晰的样式；
• 如果掩码区域有细碎小网格，可借助scipy.ndimage.morphology.binary_closing做平滑处理，过滤孤立的小区域（可选操作）。

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