做气候数据分析的人基本绕不开CMIP6这批数据。第6次耦合模式比较计划CMIP6产出的全球模拟结果在气候变化评估、区域响应、极端事件归因这些方向几乎成了标配数据源。而接触CMIP6之后你最先碰到的往往不是模式物理过程本身而是怎么用Python把NC文件打开把里面的tas变量取出来。tas也就是近地表气温2米气温是我处理得最多的变量之一。这篇就以它为例子把从读取、检查、裁剪、聚合到最终出图的完整流程走一遍。写这篇东西的受众很明确刚拿到CMIP6数据、在Jupyter Notebook里对着.nc文件不知道该点什么的同学或者已经在用但总在时间解码、经度范围、季节统计这些环节卡壳的人。文中没有高深理论绝大多数都是能直接拷贝运行的代码以及一些我踩过坑之后才总结出来的经验。CMIP6数据体量不小单文件动辄几个GB方法不对很容易卡到怀疑人生所以我把每一步为什么这么做也一并讲清楚。1. 先搞明白你要处理的东西CMIP6与tas1.1 CMIP6数据是什么文件命名怎么看CMIP6全称是第六次耦合模式比较计划参与的模式有几十个比如ACCESS-CM2、CESM2、MPI-ESM1-2-HR、CanESM5这些。每个模式都会跑一堆情景试验历史时期用historical未来预估用SSP1-2.6、SSP2-4.5、SSP5-8.5等等。这些模拟结果被打包成NetCDF格式发布也就是我们说的NC数据。下载下来之后文件名长这样tas_Amon_ACCESS-CM2_historical_r1i1p1f1_gn_185001-201412.nc这个名字不是随便起的每个字段都有含义。tas说的是变量名Amon表示大气月平均数据ACCESS-CM2是模式名historical是试验名r1i1p1f1是变体标签gn表示模式原生网格最后的185001-201412是时间范围。拿这个文件名跟别人交流时一说r1i1p1f1大家就明白你用的是第一组初始化和物理参数方案方便复现。有一个点很容易被忽略同一个模式同一套情景数据往往按时间段拆成多个文件。比如ACCESS-CM2的历史试验可能分成185001-194912、195001-199912、200001-201412三段。做长序列分析时需要先把这些文件拼起来后面我会专门讲拼接的方法。1.2 tas变量一个看起来简单但容易出错的变量tas的全称是Near-Surface Air Temperature即近地表气温通常指离地面2米处的空气温度。它几乎在所有模式的所有情景里都有输出所以特别适合作为数据集之间对比的基准变量。但你别以为简单就掉以轻心这个变量最常见的坑就是单位问题。CMIP6体系里tas的标准单位是开尔文K不是摄氏度。很多模式输出的数值常年徘徊在250到320之间第一次拿到数据的人如果忘了换算画出来的图温度高得离谱。我在实际处理中见过有人直接把开尔文当成摄氏度用了结果区域平均温度“高达”290℃还浑然不觉。换算很简单减273.15就行但这件事必须摆在数据处理流程的第一步。另外tas在不同模式下可能存在细微的物理差异比如有的模式输出的是地表气温有的模式在近地层参数化上处理不同。多模式比较时这些差异会体现在结果里。一般的使用场景下只要统一用tas、统一做单位换算差异可以接受。1.3 NC文件内部结构速览NetCDF格式可以理解为一个自带说明书的多维数据盒子文件里同时存了数据本身和描述数据的信息。打开之后你会看到四个层面的东西维度、坐标、变量、属性。维度通常就是time、lat、lon这三个但有些模式还会多出lev或plev等气压层维度。坐标是维度的具体数值比如lat从-90到90lon从0到360或从-180到180time则是一串时间戳。变量是实际存放数据的数组在CMIP6里就是tas它的形状是(time, lat, lon)。属性是写在数据旁边的元信息单位、长名称、缺省值标记都在这。理解这个结构很重要因为你后续所有操作本质上都是在跟维度、坐标、变量打交道。xarray这类的工具之所以好用核心就是把维度坐标和变量组织成一个整体让“按时间选”“按经纬度切”这些操作变成一行代码而不是像传统netCDF4库那样靠手写循环遍历下标。2. 工具链选型与环境准备2.1 为什么选xarray这条技术路线处理NC数据当然有几个选择最底层的是netCDF4库再往上是xarray。如果你只想读取一个变量的部分数据用netCDF4完全可以代码也就十几行。但如果你要做区域裁剪、时间聚合、季节平均、多文件合并这些日常操作还硬用netCDF4代码量会爆炸。我在实际项目里基本不用纯netCDF4做CMIP6分析原因有三。其一xarray的标签索引太好用了写ds.sel(time2000-01)就能直接取出某个月的数据不需要先找时间坐标的下标范围再靠数值定位。其二xarray的groupby和resample是处理时间序列的利器按年聚合、按季节平均都只需要一行代码这在气候数据分析里是刚需。其三xarray后端可以使用Dask实现懒加载和分块计算面对几个GB的单文件或者几十个文件的批量处理内存压力小很多。打个比方netCDF4像手动挡驾驶感强但每一步都要自己来xarray像自动挡你只需告诉它目的地换挡、离合的事它自己处理。多数人做研究是为了得到气候结论不是研究文件格式本身所以xarray是更合理的选择。2.2 环境搭建与安装避坑我这里假设你已经装了Python最好是3.9以上的版本。不建议直接在系统全局环境里安装这些库因为你可能同时在做爬虫、机器学习、数据处理依赖版本很容易打架。建议单独建一个conda环境conda create -n cmip python3.11 conda activate cmip然后安装核心依赖conda install -c conda-forge xarray netcdf4 cftime dask matplotlib cartopy这里特别说三个安装上容易出问题的库。cftime是处理气候历法时间的必须装否则打开某些模式的tas文件时会直接报错说time解码失败。cartopy是画地图用的依赖很多地理数据包用conda装比pip装省事得多。dask是xarray做并行和分块的后端单独处理一个文件时不是必须但面对多文件拼接时它帮大忙。如果你的网络环境访问默认源比较慢可以把conda和pip都指向国内镜像。conda可以通过在.condarc里配置channels来实现pip则用-i参数临时指定源比如清华PyPI镜像。配置方法网上一搜就有一大堆教程别在这个环节耗太多时间装好能用就行。最后建议再装一个Jupyter Notebook或VS Code写气候数据处理代码时能即时看到变量结构和数组形状调试体验比纯脚本好很多。3. 数据读取与结构探查3.1 第一步打开NC文件并看懂输出拿到文件之后第一件事是打开它看里面到底存了什么。假设文件就在当前目录下import xarray as xr ds xr.open_dataset(tas_Amon_ACCESS-CM2_historical_r1i1p1f1_gn_185001-201412.nc) print(ds)执行后会输出一大段信息类似这样xarray.Dataset Dimensions: (time: 1980, lat: 145, lon: 192) Coordinates: * time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00 * lat (lat) float64 -90.0 -88.75 -87.5 ... 88.75 90.0 * lon (lon) float64 0.0 1.875 3.75 ... 356.25 358.125 Data variables: tas (time, lat, lon) float32 ... Attributes: Conventions: CF-1.7 CMIP-6.2 title: ACCESS-CM2 output prepared for CMIP6 ...这步做完你先不急着取数据而是应该花10秒钟确认三件事。第一time维度是1980个月对应1850年1月到2014年12月正好165年。第二lat范围是-90到90lon范围是0到358.125说明这个文件的经度是0到360这种表示法后面裁剪时要格外小心。第三数据变量只有tasshape是(time: 1980, lat: 145, lon: 192)这决定了后续计算的内存需求。一个细节是xr.open_dataset默认是懒加载只读取文件元数据和变量的引用并不会真的把所有温度数组读进内存。这一点对后续处理很重要意味着你可以先探查结构再按需操作。3.2 维度、坐标和属性的检查打开数据集后用下面几个命令可以快速看维度、坐标print(ds.sizes) print(ds.coords) print(ds.tas.attrs)ds.sizes会返回一个字典告诉你每个维度的大小。ds.coords列出所有坐标变量除了time、lat、lon可能还有height或bnds这类辅助坐标。ds.tas.attrs则会返回tas的详细属性你会看到units是Klong_name是Near-Surface Air Temperature。这些信息不是摆设。我习惯在真正计算前先打印一遍属性因为不同模式、不同变量的属性细节差异很大。比如有的模式的lat坐标是降序排列从90到-90有的模式的lon是-180到180有的模式时间会带height坐标表示tas是2米处温度。你不看属性直接算可能算出完全错误的结果还不自知。还要检查一下缺测值print(ds.tas.isnull().sum().values)CMIP6绝大多数数据不含缺测但如果你处理的是插值后的数据或某些区域数据集缺测是常客。isnull().sum()这个方法会统计NaN的个数如果数量很大后续计算前就要决定是填充还是忽略。3.3 提取tas索引、切片与缺测值检查从DataSet里取出变量最简单的方式是ds[tas]或ds.tas。得到的是一个DataArray它继承了原来的维度、坐标和属性。取某个月或某个区域的数据用sel和isel# 按坐标值选择适合“人话”条件 tas_jan_2000 ds.tas.sel(time2000-01) # 按位置选择适合批量循环遍历 tas_first ds.tas.isel(time0) # 直接切片取某段纬度和经度 tas_slice ds.tas.isel(latslice(10, 20), lonslice(20, 30))这里推荐优先使用isel做粗筛因为它的判断逻辑是纯位置计算不涉及坐标值的字符串匹配速度更快。而sel更适合时间点这种语义明确的筛选比如你想取出某一年1月的数据sel(time2000-01)非常直观。缺测值方面如果前面检查发现NaN数量不高可以直接用tas_filled ds.tas.fillna(ds.tas.mean())填充。但我不建议轻易这么做缺测的处理策略跟分析目标强相关比如做趋势分析时空缺值会造成伪信号最好还是保留NaN让后续统计函数自动跳过。4. 实操过程从原始数据到可用结果4.1 单位换算开尔文转摄氏度处理气候数据时单位换算是流程里的第一步原因很简单所有后续统计、对比、可视化都依赖这个基准。CMIP6的tas默认是开尔文转换成摄氏度就是数值减273.15tas_degc ds.tas - 273.15这里有一个容易忽略的细节做减法之后新DataArray的units属性还是原来的K因为属性不会自动更新。画图时colorbar上的单位会写错或者被别人拿到数据时产生误解。所以换完单位必须同步改属性tas_degc.attrs[units] degC tas_degc.attrs[long_name] Near-Surface Air Temperature实操里还有一种情况某些数据源给的已经是摄氏度你减273.15反而错。所以每次拿到tas第一件事就是看ds.tas.attrs[units]根据实际单位决定是否转换。判断逻辑不复杂写一个小分支就行if ds.tas.attrs.get(units) K: tas ds.tas - 273.15 tas.attrs[units] degC elif ds.tas.attrs.get(units) degC: tas ds.tas else: raise ValueError(未知单位 str(ds.tas.attrs.get(units)))这样即使你批量处理很多模式的数据也能保证单位统一。4.2 区域裁剪经度范围与纬度顺序做区域研究时一般不需要全球数据裁剪出目标区域能大幅减少内存和计算量。以中国及周边区域为例通常取经度70°E到140°E、纬度15°N到55°N。但裁剪前必须先解决经度表示法的问题。我国和多数东亚研究涉及的经度是东经正数、西经负数即-180到180的表示法。而前面我们看到ACCESS-CM2的lon是0到360这种情况下直接sel(lonslice(70, 140))会取出东经70度到140度看起来没毛病但实际上遇到跨东西经的区域时就会出问题。更通用的做法是先判断是否需要把经度转换到-180到180if (ds.lon 180).any(): ds ds.assign_coords(lon(((ds.lon 180) % 360) - 180)) ds ds.sortby(ds.lon)这个公式的思路是把0到360的经度线性映射到-180到180然后按新经度排序。理解这个映射很重要经度180度本质上是同一条线0和360也是同一条线所以映射后数据本身没变只是坐标标签变了。sortby这一步必须做因为映射后坐标顺序是乱序的后续的切片和画图都依赖坐标单调递增。处理完经度再处理纬度。有的模式lat坐标是降序从90到-90如果直接sel(latslice(15, 55))因为slice要求start小于stop会取到空数据。稳妥做法是ds ds.sortby(lat)排序之后不管原始是升序还是降序都能统一处理。接下来裁剪区域tas_cn ds.tas.sel(latslice(15, 55), lonslice(70, 140))这一步之后你已经得到一份只包含目标区域的温度数据后续所有统计都在这份数据上做速度和内存都友好很多。4.3 时间聚合年际、季节与特定时期统计时间维度是CMIP6数据里最复杂的一维因为不同模式使用的日历可能不同。有的是标准公历有的是365天无闰年有的是360天日历每年固定12个月、每月30天。xarray通过cftime库自动识别日历所以时间切片基本不用操心这些差异。按月数据做成年平均用resampletas_annual tas_cn.resample(time1YS).mean(time)这里的1YS表示按年起始分组取每年所有月份的平均值。如果你要的是特定时期比如1981到2010年的气候态先切片再平均tas_recent tas_cn.sel(timeslice(1981-01, 2010-12)).mean(time)季节平均是气候分析里的高频操作用groupby可以一行搞定tas_seasonal tas_cn.groupby(time.season).mean(time)这个操作会把时间维度按MAM、JJA、SON、DJF分成四组然后每组求平均。注意groupby输出的坐标顺序不是季节的自然顺序而是MAM、JJA、SON、DJF这种字母序。如果你要单独取某季比如夏季用tas_jja tas_cn.groupby(time.season).mean(time).sel(seasonJJA)这里有个隐藏问题groupby的season是根据月份机械划分的而气候学上DJF季节横跨两个年份比如2010年12月、2011年1月和2月应该算同一个冬季。xarray的time.season并不会处理这个跨年分组所以做冬季平均时1月和2月会被归到所在年份12月会被归到下一年的DJF组里跟气候学习惯有偏差。如果你的分析对冬季定义敏感需要自己写跨年分组逻辑。好在多数区域气候分析用逐月数据直接按季节分组影响不大但要清楚这个限制。距离平均也是一个常用操作。比如分析某区域温度相对于多年平均的异常可以直接对气候态做差tas_clim tas_cn.mean(time) tas_anom tas_cn.groupby(time.year) - tas_clim这里groupby(time.year)是为了让月数据按年份分组然后整体减去气候态得到的tas_anom就是每年逐月的温度距平场。这个操作在极端事件分析里特别常用。4.4 空间平均面积权重不能省处理区域平均温度时有一个很多人忽略但影响很大的细节面积权重。地球表面的网格并不是等面积的。纬度越高网格面积越小。如果直接对lat和lon两个维度做简单平均高纬度区域的贡献会被高估热带地区的贡献会被低估。做全球平均或大区域平均时这种偏差可以达到几度足以让研究结论失真。正确的做法是按纬度余弦权重做面积加权平均。纬度余弦权重是因为球面上网格面积正比于cos(lat)weights np.cos(np.deg2rad(tas_cn.lat)) tas_region_mean (tas_cn * weights).sum(dim(lat, lon)) / weights.sum(dimlat)这个式子的含义是每个网格点的温度乘以它对应的权重然后对空间所有网格求和再除以权重总和从而得到整个区域的平均温度。注意weights是沿lat的一维数组和tas_cn的(time, lat, lon)三维数组相乘时会自动沿着lat维度广播计算逻辑正确。xarray从0.15版本开始提供了更简洁的加权平均APIregion_mean tas_cn.weighted(weights).mean(dim(lat, lon))两种写法的计算结果一样用哪种看个人习惯。我一般用第一种因为兼容性更好有时候部门服务器上的xarray版本比较老不支持weighted方法。4.5 结果导出与压缩处理完的中间结果建议保存成NC文件方便后续绘图和其他分析不用每次重头读原始大文件。保存也有讲究直接to_netcdf虽然能存但文件可能非常大因为默认不压缩。推荐加上encoding参数开启压缩tas_cn.to_netcdf( tas_cn_ACCESS-CM2_historical_1850-2014.nc, encoding{tas: {zlib: True, complevel: 4, dtype: float32}} )这里的zlib开启压缩complevel设置压缩级别数值越高文件越小但写入越慢。实际测试下来4到6这个区间性价比最高文件能缩到原来的三分之一到五分之一。dtype设为float32是因为温度精度不需要float64能省一半空间。CMIP6原始数据的精度通常就是float32保持这个精度完全没有信息损失。保存时还可以只存你真正需要的变量和坐标避免把一堆辅助坐标和属性都写进去ds_out tas_cn.to_dataset(nametas) ds_out.to_netcdf(tas_cn_ACCESS-CM2_historical_1850-2014.nc, encoding{tas: {zlib: True, complevel: 4}})这样生成的NC文件干净利落别人拿过去也能很快看懂。4.6 快速出图把处理结果画出来处理数据的最终目的通常是出图。用xarray自带的plot方法几行代码就能出一张像样的地图import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig, ax plt.subplots( figsize(8, 6), subplot_kw{projection: ccrs.PlateCarree()} ) tas_annual.sel(time2000).plot( axax, transformccrs.PlateCarree(), cmapcoolwarm, cbar_kwargs{label: °C} ) ax.coastlines(linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:) ax.set_title(Annual mean tas 2000 (ACCESS-CM2 historical)) plt.show()这里有几个容易踩的坑。第一transformccrs.PlateCarree()不能省它告诉cartopy我们的数据网格是经纬度等距网格。第二subplot_kw{projection: ccrs.PlateCarree()}指定画布的地图投影如果不指定后面的coastlines和borders都会失效。第三如果数据坐标在转换经度后映射到-180到180但画图时选择的范围不对地图上可能只显示一半这时候用ax.set_extent([70, 140, 15, 55])限制显示范围。如果只想画时间序列曲线比如区域平均温度的年际变化就更简单fig, ax plt.subplots(figsize(10, 4)) region_mean.plot(axax, colortab:red, linewidth1.2) ax.set_ylabel(Temperature anomaly (°C)) ax.set_title(Area-weighted mean tas in China region) plt.show()5. 常见问题与排查技巧实录5.1 报错速查表下面是我这几年处理CMIP6数据时遇到的典型问题整理成速查表遇到同类问题直接对号入座。报错或现象常见原因解决办法ModuleNotFoundError: No module named xarray当前Python环境没安装xarray用conda激活对应环境后安装xarray不要直接在系统环境里装KeyError: tas变量名不是tas可能是别的名称用ds.data_vars打印所有变量名确认ValueError: could not decode time缺少cftime库或日历类型不兼容安装cftime必要时指定use_cftimeTrue重新打开打开文件卡死或内存溢出数据太大懒加载被触发成了全量读入打开时加上chunks{time: 100}用Dask分块读取地图上中国区域是空白经度范围是0到360画图适配的是-180到180用assign_coords转换经度再sortby区域平均温度离谱地高或低没有做单位换算或加权平均写错检查units确认正确减273.15并用cos纬度权重多文件拼接报重复时间坐标时间段边界重叠用open_mfdataset(..., combinenested, concat_dimtime)处理画图时横坐标时间刻度叠成一团时间序列太长matplotlib自动刻度过密手动设置主刻度间隔每10年或20年显示一个标签这些坑里最常被问到的就是时间解码和经度范围后面的小节再展开。5.2 几个让我印象深刻的坑第一个是360天日历的问题。有些模式使用360天日历一年12个月每个月正好30天时间坐标上不会有1月31日。xarray处理这类时间轴时底层用cftime对象很多笨办法会在这里卡住。如果你把时间坐标直接转成pandas的DatetimeIndex九成会报错。正确做法是保持cftime类型别乱转。做时间分析时resample和groupby在xarray里对cftime是兼容的所以尽量别自己手动操作时间索引。第二个是纬度顺序颠倒。某个模式的lat坐标是降序从90到-90我没检查直接按区域裁剪sel(latslice(15, 55))返回空数组当时愣了好一会儿才反应过来。从那以后我只要打开新文件就会先打印ds.lat.values[:5]和ds.lat.values[-5:]确认顺序。这也是为什么我在前面代码里强调sortby(lat)这步不能省。第三个是多文件拼接时的时间边界重叠。有些ESGF节点下载的数据时间段是重合的比如一个文件到2005年12月另一个文件从2005年1月开始open_mfdataset默认按坐标合并就会出现重复坐标报错。解决办法是显式指定合并方式ds xr.open_mfdataset( tas_Amon_ACCESS-CM2_historical_r1i1p1f1_gn_*.nc, combinenested, concat_dimtime, data_varsminimal, coordsminimal, parallelTrue )combinenested表示把文件按顺序拼接而不是按坐标自动取交集data_varsminimal和coordsminimal避免拼接时重复处理非索引坐标。这样即使时间边界有一点重叠也能拼成功。第四个让我印象很深的坑跟面积权重有关。早期处理全球平均温度时我直接tas.mean(dim(lat, lon))结果跟官方发布的全球平均温度序列对不上差了不少。后来翻资料才发现CMIP6分析社区做全球平均时普遍使用面积权重而且有些研究还强调要用更精确的quadrature权重。从那以后我养成了一个习惯不管做区域平均还是全球平均先问自己一句“面积权重加了没”。也建议你把这个习惯刻在脑子里。还有一个跟数据量有关的经验。CMIP6单文件经常是1到3GB处理时临时文件、中间结果会很快把磁盘占满。我自己的习惯是下载和处理的原始文件放在单独的数据盘处理完的中间产物及时清理最后只保留压缩后的结果文件。毕竟一个模式一个情景就有好几GB几十个模式处理下来数据管理比代码本身更考验人。最后分享一个小技巧如果你要处理的文件时间跨度特别长比如1850到2100年可以先做裁剪再聚合不要上来就全量载入内存。配合chunks参数和Dask的懒加载哪怕是几十GB的数据也能在普通办公电脑上跑完。这也是xarray这套工具链真正强大的地方。