我从ESGF节点折腾了一整天下载完某个模式的历史模拟数据满心欢喜地用Python打开NetCDF文件结果时间轴是noleap日历经度坐标居然不是规则网格降水变量单位是kg m-2 s-1而不是我熟悉的mm/day。这种滋味做过CMIP6数据处理的朋友应该都不陌生。CMIP6第六次国际耦合模式比较计划是当前气候变化研究最核心的数据源全球几十个模式组贡献了海量模拟结果支撑了IPCC第六次评估报告。但很多人第一次拿到CMIP6数据时都会被它的复杂性劝退——文件命名规则、实验场景体系、网格结构、变量单位、偏差特征每一个环节都有坑。本文不打算重复官方文档而是从我个人的科研实操经验出发完整梳理CMIP6数据从下载、预处理到在气候变化、水文、生态领域实际落地的全过程把那些文档里不会写的细节一并抖出来。1. CMIP6数据到底特殊在哪规模、结构与命名规则很多新手第一次打开CMIP6数据时脑子里还是CMIP5时代的习惯。但实际上CMIP6在数据组织和复杂度上做了相当大的改动如果不先把这些底层逻辑搞清楚后面每一步都会出问题。1.1 数据规模和文件体系几十TB起步的资源池CMIP6总数据量比CMIP5提升了数倍。全球各模式中心产出的数据通过ESGFEarth System Grid Federation地球系统网格联盟分布式节点发布。单个变量的单次实验数据动辄几十GB到数百GB比如pr降水日数据一个模式、一个情景、一个成员全球范围0.5°~2°分辨率时间跨度几十年压缩后的NetCDF文件也有数十GB。这个规模直接决定了处理策略。我的习惯是坚决避免把整个文件load()进内存全程使用xarraydask的分块懒加载模式。后面第三章会细说。1.2 文件命名规则一串字符把家底全交代了CMIP6的标准文件名长这样pr_day_ACCESS-CM2_historical_r1i1p1f1_gn_20100101-20141231.nc逐个拆解字段含义典型取值pr变量名tas, pr, tasmax, tasmin, sfcWind…day时间频率day, mon, 3hr, 6hr…ACCESS-CM2模式名称全球约100多个模式historical实验名称historical, ssp245, ssp585…r1i1p1f1成员编号初始化方案r、扰动参数i、物理方案p、强迫方案fgn网格类型gn原生网格, gr规则网格20100101-20141231时间覆盖范围每个文件通常覆盖5-10年这个命名规则里有几个关键点。第一historical实验一般覆盖1850年到2014年ssp245、ssp585等情景实验从2015年开始两者是连续的做长期趋势分析时要拼接。第二r1i1p1f1不同成员之间存在气候内部变率差异如果你只下一个成员得到的结果可能带有特定的内部变率信号。对于某些研究比如极端事件频率分析建议至少取3-5个成员做集合。第三gn原生网格意味着经纬度坐标不是规则递增的不同模式可能有不同的网格结构高斯网格、立方球网格等后续插值必须做。还有一点容易被忽略不同模式的历史模拟时段并不完全一致有的从1850年开始有的从1851年开始。合并多个模式时一定先取公共时间范围。1.3 实验场景与情景路径不是只有SSP1-2.6和SSP5-8.5CMIP6有一个庞大的情景矩阵ScenarioMIP核心是共享社会经济路径SSP与代表性浓度路径RCP的组合。最常用的四个SSP1-2.6可持续发展低温室气体排放2100年辐射强迫约2.6 W/m²SSP2-4.5中间路径约4.5 W/m²SSP3-7.0区域竞争高排放约7.0 W/m²SSP5-8.5化石燃料密集型发展约8.5 W/m²此外还有ssp126、ssp370、ssp119等变体。水文和生态研究用到ssp126和ssp370的频率也很高。我个人的建议是做研究时不要只看SSP5-8.5这一个高排放情景多情景对比才能真正反映未来不确定性范围。IPCC第六次评估报告的核心结论很多是基于中等情景SSP2-4.5得出的。2. 从下载到归档ESGF使用体验与本地数据规划获取CMIP6数据看似简单——打开ESGF节点搜索下载就行。但实际操作中下载策略不当会浪费大量时间。2.1 ESGF节点选择与下载方式别用浏览器手动点ESGF是全球分布式网络全球有多个镜像节点。国内用户可能对国外节点的访问速度有体会全文不展开网络工具的讨论。我推荐的做法是直接用wget脚本批量下载而不是在网页上一个一个点。ESGF搜索界面返回结果后可以选择List All然后生成wget脚本服务器会返回一个.sh文件里面有所有文件的下载链接和MD5校验信息。# 下载脚本的基本用法 bash download_script.sh脚本默认只下载当前节点上存在的文件如果某个文件在节点间同步失败它会自动跳过。还有几个参数值得注意--checksum强制做MD5校验建议加上防止数据损坏--timeout设置超时时间默认进程数可能只有3可以用环境变量调大下载后先做一次文件完整性检查md5sum -c checksum_file.txt integrity_check.log2.2 本地目录结构与元数据管理CMIP6数据量大、文件多没有良好的目录规划找到某个具体文件都要花半天。我推荐按以下结构归档CMIP6/ └── CMIP/ └── MRI-ESM2-0/ ├── historical/ │ └── r1i1p1f1/ │ ├── day/ │ │ └── pr/ │ └── mon/ │ └── tas/ └── ssp245/ └── r1i1p1f1/ ├── day/ └── mon/这种结构跟ESGF节点的目录结构保持一致后续要批量遍历文件时非常省事。同时我建议在本地维护一个metadata.csv记录每个文件对应的模式、实验、变量、频率、成员、时间范围和文件路径。用Python的glob扫一遍目录加上文件名解析几百行代码就能搞定一个简易的本地数据目录索引。import glob import pandas as pd files glob.glob(/data/CMIP6/**/*.nc, recursiveTrue) records [] for f in files: fname f.split(/)[-1] parts fname.replace(.nc, ).split(_) if len(parts) 7: variable, freq, model, experiment, member, grid, time_range parts[:7] records.append({ variable: variable, freq: freq, model: model, experiment: experiment, member: member, grid: grid, time_range: time_range, file_path: f, }) df pd.DataFrame(records) df.to_csv(metadata.csv, indexFalse)这个索引文件后续做多模式集合、批量处理时能节省大量时间。3. 用xarray重构CMIP6预处理流水线CMIP6数据处理通常涉及格式检查、时间轴处理、空间裁剪、网格插值、单位换算等多个步骤。过去很多人习惯用NCO或CDO命令行工具但对于科研人员来说用Python的xarray生态做全流程处理更灵活也更容易复现。这里我给出一个常用流水线的核心代码框架。3.1 读取与时间轴修复cftime日历是第一道坎打开CMIP6文件的第一步就必须处理日历问题。CMIP6的NetCDF文件遵循CF约定常见日历有noleap一年恒为365天、365_day、360_day等。直接用pandas处理会直接报错或产生错位。import xarray as xr import cftime ds xr.open_dataset( pr_day_ACCESS-CM2_historical_r1i1p1f1_gn_20100101-20141231.nc, use_cftimeTrue ) # 查看时间维信息 print(ds.time)关键是把时间轴转换为标准的公历日历。xarray中可以用cftime提供的date2num或者直接使用ds.time的dt接口。常见的做法是转为datetime64# 如果是noleap日历直接转换为标准日历 ds[time] ds.indexes[time].to_datetimeindex()to_datetimeindex()在数据量大时可能较慢但它能把noleap的365天虚拟日映射到标准日期上。需要注意的是noleap日历的12月31日并不存在转换后时间轴末尾可能略有偏差在做月统计时务必检查最后几个时间点。读取后建议先做一个数据体检检查变量维度、坐标、缺失值比例def quick_inspect(ds): print(fData variables: {list(ds.data_vars)}) for var in ds.data_vars: da ds[var] print(f\n{var}: shape{da.shape}, dtype{da.dtype}) print(f attrs: units{da.attrs.get(units, N/A)}) # 简单缺失值检查 if da.dtype.kind in fc: nan_pct (da.isnull().sum() / da.size).values print(f NaN ratio: {nan_pct:.2%})这一步能帮你提前发现问题避免后续计算到一半才发现数据有问题。3.2 单位统一与变量换算CMIP6的变量遵循严格的标准名和单位。常见变量变量名含义CMIP6原始单位常用工程单位tas / tasmax / tasmin近地面气温K°Cpr降水含液态固态kg m⁻² s⁻¹mm/daysfcWind近地面风速m/sm/srsds地表下行短波辐射W/m²MJ/m²/daymrro地表径流含baseflow或总径流kg m⁻² s⁻¹mm/monthevspsbl蒸散发kg m⁻² s⁻¹mm/day温度换算是最简单的减273.15即可。降水要乘以86400才能从kg m⁻² s⁻¹换算成mm/day因为1 kg/m²的水深正好是1 mm而1天有86400秒。很多新手忘记这一步直接把模式降水跟站点观测对比结果数值小到离谱。ds[tas] ds[tas] - 273.15 ds[tas].attrs[units] degC ds[pr_mmday] ds[pr] * 86400 ds[pr_mmday].attrs[units] mm/day3.3 时间重采样与空间裁剪做月尺度分析时日数据需要聚合成月数据。xarray的resample功能非常方便但要注意聚合时对time_bnds的处理——聚合后的值应该是加权平均还是简单平均取决于变量性质。对降水来说月总量应该是日值的和对温度来说月均应该是日值的平均。# 月降水总量mm/month pr_monthly ds[pr_mmday].resample(time1M).sum() # 月均温度°C tas_monthly ds[tas].resample(time1M).mean()空间裁剪也有两种常见方式。一种是简单的经纬度边界裁剪ds_region ds.sel(lonslice(80, 130), latslice(3, 53))这种方式适合研究区域是矩形的情况。如果研究区是流域或行政边界推荐用regionmask包做掩膜import regionmask # 以中国省级边界为例 regions regionmask.defined_regions.natural_earth_v5_0_0.countries_110 mask regions.mask(ds_region[lon], ds_region[lat]) ds_masked ds_region.where(mask 46) # 中国对应的编号3.4 网格插值多模式集合必须迈过的坎CMIP6不同模式的空间分辨率差异很大有的原生网格是2°见方有的高达0.5°。做多模式集合平均或与观测数据对比时必须统一到同一网格。xesmf是目前最顺手的重网格工具底层调用ESMF库接口简洁。以双线性插值为例import xesmf as xe import xarray as xr # 目标网格1°×1°规则网格 grid_out xr.Dataset( { lat: ([lat], np.arange(-90, 90, 1.0)), lon: ([lon], np.arange(0, 360, 1.0)), } ) # 确定源网格 grid_in ds_regridded.copy() # 构建回归器双线性插值 regridder xe.Regridder(grid_in, grid_out, bilinear) # 执行插值 ds_out regridder(ds_regridded, keep_attrsTrue)注意几点第一CMIP6的gn原生网格数据必须先转换到规则网格。实际上ESGF网站上也可以选择下载gr规则网格版本如果不需要处理原生网格特性直接下载gr版本能省掉很多事情。第二插值本身会引入平滑效应对于降水这种空间变异性强的变量插值后的极值是有偏差的。如果做极端降水分析我更推荐保守插值conservative能更好保持区域总量不过计算量会大一些。第三插值权重只与网格有关可以保存下来重复使用regridder.to_netcdf(ACCESS_1x1_bilinear_weights.nc) # 下次直接加载 regridder xe.Regridder(None, None, reuse_weightsTrue)4. 气候变化、水文、生态三个领域的落地打法数据预处理完不同领域的应用路数差异很大。这里分别给出最常用的分析框架。4.1 气候变化研究趋势、极端指数与不确定性量化最基础的工作是分析气温和降水的变化趋势、突变点及极端事件频率变化。经典做法包括Mann-Kendall趋势检验和Sens slope估计。scipy配合pymannkendall包即可完成。import pymannkendall as mk import numpy as np # 假设annual_temp是年序列 trend, h, p, z, Tau, s, var_s, slope, intercept mk.original_test(annual_temp)极端气候指数的计算推荐直接使用xclim库它已经把ETCCDI的全部指数都实现好了。例如计算夏季日最高温超过某阈值的天数TX90p或连续干旱日数CDDimport xclim.indices as xci # CDD日降水1mm的连续最长天数 cdd xci.max_consecutive_dry_days(ds[pr_mmday], thresh1.0) # TXx年最大日最高温 txx xci.tx_max(ds[tasmax])多模式集合分析时一定要区分两种不确定性来源模式间差异不同模式对同一物理过程的参数化不同和内部变率同一模式不同初始条件导致的结果差异。一般用多模式集合平均MME降低内部变率影响用模式间标准差spread表征不确定性范围。4.2 水文领域把CMIP6输出喂给水文模型水文模拟通常需要降尺度和偏差校正后的气象驱动数据。这里有几个常用的处理路径。偏差校正是水文应用的关键步骤因为GCM输出的降水频率和强度分布与实测有明显系统偏差。最简单的Delta方法适合温度obs_clim obs.groupby(time.month).mean() # 观测月气候态 gcm_clim gcm.groupby(time.month).mean() # 模式历史月气候态 # 温度加性校正 tas_corrected tas - gcm_clim obs_clim降水的Delta方法用乘法比率但在月降水接近0时容易产生异常值更推荐分位数映射Quantile Mapping, QM。scipy的quantile或直接按经验分布函数做映射即可实现这里给个简化版from scipy.stats import norm def quantile_mapping(obs, model_hist, model_future): 分位数映射纠偏将model_future的分布映射到obs的分布 n len(obs) obs_sorted np.sort(obs) model_hist_sorted np.sort(model_hist) # 经验CDF插值 model_fut_uniform np.interp(model_future, model_hist_sorted, np.linspace(0, 1, n)) corrected np.interp(model_fut_uniform, np.linspace(0, 1, n), obs_sorted) return corrected实际项目里建议用python-cmethods或xclim.sdba支持多种成熟的偏差校正算法包括分位数增量映射、日际循环重构等。在水文领域模式数据一般还需要降尺度到水文模型网格如5 km或1 km降尺度方法的选择会显著影响径流模拟结果。4.3 生态领域物候、生产力与碳循环模拟生态应用更关注温度、降水与生物过程的关系讲究的是驱动数据的连续性和协调性。物候分析通常需要连续无缺测的温度序列。基于GDD生长度日的物候模型直接采用日平均温度模式数据在此问题上表现相对稳定。不过要注意模式温度的长期漂移可能导致假物候变化我建议先做整体偏差校正再驱动物候模型。植被生产力模拟时如果直接使用GCM输出的温度和降水驱动生态过程模型如Biome-BGC、LPJ-GUESS原始数据的时间分辨率、辐射变量和湿度变量的完整性非常关键。CMIP6里这些变量分布在不同MIP中例如LUMIP提供土地利用变化数据C4MIP提供碳循环反馈数据。做碳循环相关研究时需要检查模式是否耦合了动态植被和碳循环模块不同模式对CO₂施肥效应的模拟差异很大。一个容易忽略的细节是CMIP6模式模拟的是潜在植被还是实际土地利用下的状态不同模式定义不同生态应用里要仔细读mip_era和experiment说明文档不能想当然。5. 那些年我们踩过的CMIP6数据坑CMIP6数据处理的坑多到可以单独写一本书。我把自己和身边同事踩过的典型坑整理一下给大家提个醒。5.1 日历与时间戳的隐性错误noleap日历的天数编码方式与标准日历不同直接用datetime64[ns]读取时可能整体偏移。我遇到过一次情况某模式数据的时间戳是1960-01-01开始编码的days since格式转成公历时所有日期都错位了一天导致月平均结果全部错误。排查方法是单独提取时间轴做日和月两个尺度的交叉检验。另一个常见问题是time_bnds缺失。某些模式文件的time_bnds含日界信息但个别数据版本没有聚合计算年值时会引入半个格点的窗口误差。建议无论原始文件是否带time_bnds都在处理第一环节用ds.time.to_index().to_period(D)重新构建边界并检查。5.2 变量单位与因子换算的灵魂拷问除了温度和降水其他变量也有各自的换算陷阱。比如evspsbl蒸散发的单位是kg m-2 s-1换算成mm/day同样是乘以86400。但注意CMIP6中部分模式把evspsbl定义为参考蒸散发而不是实际蒸散发两者含义完全不同。拿到数据先看cell_methods属性里的说明再决定如何参与计算。rsds和rsus上行/下行短波辐射的符号约定、净辐射的正负号在不同模式下有过不统一的历史问题CMIP6已经统一了方向但仍建议每个模式都检查一遍极值范围是否合理。5.3 实验成员乱选导致的结果偏差很多人在ESGF下载时随意选了一个r1i1p1f1觉得都是同一个模式。实际上不同初始条件会带来很大的内部变率差异尤其对于区域尺度的极端事件这种差异可能接近甚至超过模式间的差异。个人经验是做区域气候风险评估至少用5个成员的集合平均否则不确定性的估计是相当不可靠的。好在CMIP6大部分模式的r1i1p1f1~r5i1p1f1的数据都齐全批量下载处理也可以流水线化。5.4 处理效率瓶颈一个一个文件处理真的是给自己挖坑如果用xarray直接读取几十GB的大文件不做分块优化内存很容易爆。正确做法是设置chunks参数让dask把计算拆成小块并行处理ds xr.open_dataset( pr_day_ACCESS-CM2_historical_r1i1p1f1_gn_20100101-20141231.nc, chunks{time: 3650, lat: 100, lon: 100} )同时推荐把多个时间切片文件先合并成单一大文件xr.open_mfdataset做一次NetCDF4的压缩存储后续读取效率会高很多。因为每次开文件都有I/O开销文件数量越多重复读取的损耗越大。# 合并多时间切片 ds_merged xr.open_mfdataset( /data/CMIP6/CMIP/*/historical/r1i1p1f1/day/pr/*.nc, combineby_coords, chunksauto )合并之后可以输出成一个年份段更长的文件或者转成Zarr格式。Zarr在读写效率上比NetCDF更好后续增量计算很方便。对于要做大量重采样和集合操作的场景这个改动带来的速度提升非常明显。写在最后CMIP6数据是个庞然大物但也是气候变化、水文和生态研究的富矿。从我个人的体验来说做这套数据处理最核心的其实是三件事把命名规则搞透彻把时间日历和单位换算这类基础但致命的细节固定成标准流程以及建立一套高效可复现的代码框架。很多人研究做不下去其实不是分析思路的问题而是数据处理环节就已经耗尽了精力。最后分享一个小技巧不要每次拿到新数据就从零写处理代码建议把数据体检、单位换算、时间重采样、网格插值封装成几个标准函数用yaml配置文件把模式和实验参数固化下来形成自己团队的CMIP6专属流水线。这样无论是换一个模式、还是加一个情景研究改两行配置就能跑通可以节约大量时间也大幅降低出错概率。
