<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排序,便于比较模式性能