<think>
我们需要处理24个CMIP6模式，这些模式被分为8类，每类有不同的路径和命名规则。
步骤：
1. 定义每个模式的信息：包括模式名称、路径、文件名格式（注意不同类别的文件名中实验ID和网格标签不同）
2. 循环1980-2014年，对每个模式：
a. 读取观测数据（每年一个文件）
b. 读取模式数据（每年一个文件）
c. 将模式数据从K转换为℃（减去273.15）
d. 裁剪到中国区域（18-54N, 73-135E）
e. 计算年平均（时间维平均）
f. 计算该模式该年在中国区域的RMSE（空间上每个格点先计算年平均值，然后计算空间上所有格点的RMSE）
3. 对每个模式，计算1980-2014年RMSE的平均值（时间平均）
4. 绘制柱状图，每个模式一个柱子，柱子高度为多年平均RMSE
注意：由于有24个模式，我们需要为每个模式存储一个列表，记录每年的RMSE，然后取平均。但我们最终展示的是每个模式在1980-2014年的平均RMSE。
由于模式众多，我们首先构建一个模式信息列表，每个元素是一个字典，包含：
- 模式名称（用作标识）
- 路径
- 文件名模板（用模式名称替换XX，年份替换xxxx）
- 注意：不同类别的模板中，实验ID（如r1i1p1f1）和网格标签（gn/gr/gr1）不同，所以需要分别指定。
根据描述，将24个模式按照类别分为：
第一类（10个）：ACCESS-CM2, ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, MIROC6, MPI-ESM1-2-HR, MPI-ESM1-2-LR, MRI-ESM2-0, CMCC-ESM2
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/XX
文件名模板：tas_day_{model}_historical_r1i1p1f1_gn_{year}.nc
第二类（2个）：CESM2-WACCM, FGOALS-g3
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/YY
文件名模板：tas_day_{model}_historical_r3i1p1f1_gn_{year}.nc
第三类（3个）：EC-Earth3, EC-Earth3-Veg-LR, IPSL-CM6A-LR
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/ZZ
文件名模板：tas_day_{model}_historical_r1i1p1f1_gr_{year}.nc
第四类（5个）：GFDL-CM4, GFDL-ESM4, INM-CM4-8, INM-CM5-0, KIOST-ESM
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/AA
文件名模板：tas_day_{model}_historical_r1i1p1f1_gr1_{year}.nc
第五类（2个）：CNRM-CM6-1, CNRM-ESM2-1
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/BB
文件名模板：tas_day_{model}_historical_r1i1p1f2_gr_{year}.nc
第六类（1个）：CESM2（注意第二类中已经有CESM2-WACCM，这里是CESM2）
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/CESM2
文件名模板：tas_day_CESM2_historical_r4i1p1f1_gn_{year}.nc
注意：这里模式名称是固定的，所以我们可以用单独的一个条目。
第七类（1个）：MIROC-ES2L
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/MIROC-ES2L
文件名模板：tas_day_MIROC-ES2L_historical_r1i1p1f2_gn_{year}.nc
第八类（1个）：GISS-E2-1-G
路径：/media/HUAWEI/Elements SE/ec-adjusted/historical/GISS-E2-1-G
文件名模板：tas_day_GISS-E2-1-G_historical_r1i1p1f2_gn_{year}.nc
注意：在构建模式列表时，我们为每个模式指定其路径和文件名格式。
由于每个模式的文件命名规则不同，我们按照上述分类构建模式信息列表。
另外，观测数据的路径和命名是固定的。
我们将为每个模式计算1980-2014年每年的RMSE，然后对年份求平均，得到该模式的平均RMSE。
最后，绘制柱状图，横轴为模式名称，纵轴为平均RMSE。
由于模式名称较长，可能需要旋转横轴标签。
代码结构：
1. 定义常量（观测路径、输出路径、中国区域等）
2. 构建模式信息列表（每个元素包含：模式名，路径，文件名模板）
3. 循环每个模式
初始化一个列表存储该模式每年的RMSE（35年）
循环1980-2014年：
读取观测数据（如果某年观测数据不存在，则跳过该年；同样，模式数据不存在也跳过）
读取模式数据（注意文件名模板中的年份和模式名）
处理数据（单位转换，区域选择，计算年平均）
计算该年的RMSE（空间上）
记录该年的RMSE
计算该模式多年平均RMSE（跳过缺失年份）
4. 收集每个模式的平均RMSE
5. 绘制柱状图
注意：由于数据量较大，我们使用xarray的延迟加载。但是循环35年*24个模式，可能会占用较多内存，注意及时关闭文件并清理内存。
我们将使用上下文管理器打开文件，确保文件被关闭。
考虑到可能有些年份缺失，我们使用try-except跳过。
由于模式数量多，我们打印进度信息。
实现：
</think>以下是为24个CMIP6模式计算1980-2014年中国区域气温模拟RMSE并绘制柱状图的完整代码：
```python
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pathlib import Path
import os
# 设置路径
OBS_PATH = Path("/media/HUAWEI/LY-TSD302/ChinaMet_0.25deg(365_fixed)")
OUTPUT_DIR = Path("/media/HUAWEI/Elements SE/1111")
os.makedirs(OUTPUT_DIR, exist_ok=True)
# 中国区域范围
CN_REGION = {"lon": slice(73, 135), "lat": slice(18, 54)}
# 定义所有CMIP6模式的信息 (名称, 路径, 文件名格式)
MODELS = [
# 第一类: r1i1p1f1_gn
{"name": "ACCESS-CM2", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/ACCESS-CM2",
"pattern": "tas_day_ACCESS-CM2_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "ACCESS-ESM1-5", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/ACCESS-ESM1-5",
"pattern": "tas_day_ACCESS-ESM1-5_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "BCC-CSM2-MR", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/BCC-CSM2-MR",
"pattern": "tas_day_BCC-CSM2-MR_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "CanESM5", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CanESM5",
"pattern": "tas_day_CanESM5_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "MIROC6", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/MIROC6",
"pattern": "tas_day_MIROC6_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "MPI-ESM1-2-HR", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/MPI-ESM1-2-HR",
"pattern": "tas_day_MPI-ESM1-2-HR_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "MPI-ESM1-2-LR", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/MPI-ESM1-2-LR",
"pattern": "tas_day_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "MRI-ESM2-0", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/MRI-ESM2-0",
"pattern": "tas_day_MRI-ESM2-0_historical_r1i1p1f1_gn_{year}.nc"},
{"name": "CMCC-ESM2", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CMCC-ESM2",
"pattern": "tas_day_CMCC-ESM2_historical_r1i1p1f1_gn_{year}.nc"},
# 第二类: r3i1p1f1_gn
{"name": "CESM2-WACCM", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CESM2-WACCM",
"pattern": "tas_day_CESM2-WACCM_historical_r3i1p1f1_gn_{year}.nc"},
{"name": "FGOALS-g3", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/FGOALS-g3",
"pattern": "tas_day_FGOALS-g3_historical_r3i1p1f1_gn_{year}.nc"},
# 第三类: r1i1p1f1_gr
{"name": "EC-Earth3", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/EC-Earth3",
"pattern": "tas_day_EC-Earth3_historical_r1i1p1f1_gr_{year}.nc"},
{"name": "EC-Earth3-Veg-LR", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/EC-Earth3-Veg-LR",
"pattern": "tas_day_EC-Earth3-Veg-LR_historical_r1i1p1f1_gr_{year}.nc"},
{"name": "IPSL-CM6A-LR", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/IPSL-CM6A-LR",
"pattern": "tas_day_IPSL-CM6A-LR_historical_r1i1p1f1_gr_{year}.nc"},
# 第四类: r1i1p1f1_gr1
{"name": "GFDL-CM4", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/GFDL-CM4",
"pattern": "tas_day_GFDL-CM4_historical_r1i1p1f1_gr1_{year}.nc"},
{"name": "GFDL-ESM4", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/GFDL-ESM4",
"pattern": "tas_day_GFDL-ESM4_historical_r1i1p1f1_gr1_{year}.nc"},
{"name": "INM-CM4-8", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/INM-CM4-8",
"pattern": "tas_day_INM-CM4-8_historical_r1i1p1f1_gr1_{year}.nc"},
{"name": "INM-CM5-0", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/INM-CM5-0",
"pattern": "tas_day_INM-CM5-0_historical_r1i1p1f1_gr1_{year}.nc"},
{"name": "KIOST-ESM", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/KIOST-ESM",
"pattern": "tas_day_KIOST-ESM_historical_r1i1p1f1_gr1_{year}.nc"},
# 第五类: r1i1p1f2_gr
{"name": "CNRM-CM6-1", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CNRM-CM6-1",
"pattern": "tas_day_CNRM-CM6-1_historical_r1i1p1f2_gr_{year}.nc"},
{"name": "CNRM-ESM2-1", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CNRM-ESM2-1",
"pattern": "tas_day_CNRM-ESM2-1_historical_r1i1p1f2_gr_{year}.nc"},
# 第六类: CESM2
{"name": "CESM2", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/CESM2",
"pattern": "tas_day_CESM2_historical_r4i1p1f1_gn_{year}.nc"},
# 第七类: MIROC-ES2L
{"name": "MIROC-ES2L", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/MIROC-ES2L",
"pattern": "tas_day_MIROC-ES2L_historical_r1i1p1f2_gn_{year}.nc"},
# 第八类: GISS-E2-1-G
{"name": "GISS-E2-1-G", "path": "/media/HUAWEI/Elements SE/ec-adjusted/historical/GISS-E2-1-G",
"pattern": "tas_day_GISS-E2-1-G_historical_r1i1p1f2_gn_{year}.nc"}
]
# 存储所有模式的结果
all_results = {}
# 循环处理所有模式
for model_info in MODELS:
model_name = model_info["name"]
print(f"\n开始处理模式: {model_name}")
# 存储当前模式的结果
model_results = {"year": [], "rmse": []}
years_processed = 0
# 循环处理1980-2014年
for year in range(1980, 2015):
# 观测文件
obs_file = OBS_PATH / f"ChinaMet_025deg_tmpmean_{year}.nc"
# 模式文件
model_file = Path(model_info["path"]) / model_info["pattern"].format(year=year)
# 检查文件是否存在
if not os.path.exists(obs_file):
print(f" 观测数据缺失: {obs_file}")
continue
if not os.path.exists(model_file):
print(f" 模式数据缺失: {model_file}")
continue
try:
# 读取并处理观测数据
with xr.open_dataset(obs_file) as ds_obs:
obs_annual = ds_obs['tmpmean'].sel(**CN_REGION).mean(dim='time')
# 读取并处理模式数据
with xr.open_dataset(model_file) as ds_model:
# 查找温度变量名 (通常是'tas'，但确保正确)
temp_var = 'tas'
if 'tas' not in ds_model:
# 尝试查找其他可能变量名
possible_vars = [v for v in ds_model.data_vars if v.lower().startswith('tas')]
if possible_vars:
temp_var = possible_vars[0]
else:
raise ValueError(f"在{model_file}中未找到温度变量")
model_annual = (ds_model[temp_var] - 273.15).sel(**CN_REGION).mean(dim='time')
# 计算空间均方根误差 (RMSE)
mse = np.nanmean((model_annual - obs_annual)**2)
rmse = np.sqrt(mse)
# 保存结果
model_results["year"].append(year)
model_results["rmse"].append(rmse)
years_processed += 1
print(f" {year}年 RMSE: {rmse:.4f}℃")
except Exception as e:
print(f" 处理{year}年数据出错: {str(e)}")
continue
# 计算多年平均RMSE
if years_processed > 0:
mean_rmse = np.mean(model_results["rmse"])
print(f"{model_name} 平均RMSE ({years_processed}年): {mean_rmse:.4f}℃")
# 存储到总结果
all_results[model_name] = {
"mean_rmse": mean_rmse,
"years_processed": years_processed,
"yearly_rmse": model_results["rmse"]
}
else:
print(f"{model_name} 无有效年份数据")
all_results[model_name] = {
"mean_rmse": np.nan,
"years_processed": 0,
"yearly_rmse": []
}
# 准备绘图数据
model_names = []
mean_rmses = []
years_counts = []
for model, data in all_results.items():
model_names.append(model)
mean_rmses.append(data["mean_rmse"])
years_counts.append(data["years_processed"])
# 创建DataFrame用于排序和绘图
results_df = pd.DataFrame({
"Model": model_names,
"Mean_RMSE": mean_rmses,
"Years_Count": years_counts
})
# 按RMSE从小到大排序
results_df = results_df.sort_values(by="Mean_RMSE", ascending=True)
# 绘制柱状图
plt.figure(figsize=(18, 10))
bars = plt.bar(results_df["Model"], results_df["Mean_RMSE"],
color=plt.cm.viridis(np.linspace(0, 1, len(results_df))))
# 添加数值标签
for i, bar in enumerate(bars):
height = bar.get_height()
plt.text(bar.get_x() + bar.get_width()/2., height + 0.05,
f'{height:.2f}\n({years_counts[i]})',
ha='center', va='bottom', fontsize=9)
# 设置图表属性
plt.xlabel('气候模式', fontsize=14, labelpad=15)
plt.ylabel('平均RMSE (℃)', fontsize=14, labelpad=15)
plt.title('24个CMIP6模式中国区域年平均气温模拟能力评估 (1980-2014)',
fontsize=16, pad=20)
plt.xticks(rotation=45, ha='right', fontsize=11)
plt.ylim(0, max(mean_rmses) * 1.25 if not np.all(np.isnan(mean_rmses)) else 5)
# 添加网格线和图例
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.tight_layout()
# 保存图表
output_file = OUTPUT_DIR / "CMIP6_Models_RMSE_Evaluation.png"
plt.savefig(output_file, dpi=300, bbox_inches='tight')
print(f"\n评估完成! 结果图表已保存至: {output_file}")
# 保存结果到CSV文件
csv_file = OUTPUT_DIR / "CMIP6_Models_RMSE_Results.csv"
results_df.to_csv(csv_file, index=False, encoding='utf-8-sig')
print(f"结果数据已保存至: {csv_file}")
# 显示图表
plt.show()
```
### 代码功能说明：
1. **多模式处理框架**：
- 定义了24个CMIP6模式的信息（名称、路径、文件命名规则）
- 自动处理8类不同命名规则的模式数据
- 统一计算每个模式1980-2014年的RMSE
2. **智能数据读取**：
- 自动检测温度变量名称（优先使用'tas'，自动尝试匹配其他名称）
- 处理不同模式的文件命名差异（r1i1p1f1/r3i1p1f1/r4i1p1f1等）
- 自动处理网格标签差异（gn/gr/gr1）
3. **错误处理与日志**：
- 详细记录每个模式每年的处理状态
- 自动跳过缺失文件或无效数据
- 记录每个模式的有效年份数量
4. **结果分析与可视化**：
- 计算每个模式的多年平均RMSE
- 按RMSE值从小到大排序展示
- 使用渐变色区分不同模式
- 在柱状图上显示RMSE值和有效年份数
- 保存完整结果到CSV文件
### 输出结果：
1. **图表文件**：
- `/media/HUAWEI/Elements SE/1111/CMIP6_Models_RMSE_Evaluation.png`
- 展示24个模式按RMSE排序的柱状图
- 每条柱子上方显示RMSE值和有效年份数
2. **数据文件**：
- `/media/HUAWEI/Elements SE/1111/CMIP6_Models_RMSE_Results.csv`
- 包含每个模式的名称、平均RMSE和有效年份数
- 按RMSE值从小到大排序
### 使用说明：
1. 代码自动创建输出目录
2. 自动处理所有24个CMIP6模式
3. 控制台显示详细处理进度和结果
4. 最终结果包含图表和数据文件
5. 缺失年份数据会自动跳过，不影响其他处理
### 注意事项：
1. 模式数据路径需确保正确
2. 部分模式数据可能有变量名差异，代码已做兼容处理
3. 如果某模式多年数据缺失，其RMSE将显示为NaN
4. 柱状图按RMSE排序，便于比较模式性能
博客提及了逻辑相关内容，但信息较少。逻辑在信息技术领域有重要应用，如逻辑编程等。
扫一扫
447
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
