搞气象数据分析的人多多少少都体会过那种“手里攒了一堆数据却不知道怎么顺畅下手”的憋屈感。我最早处理再分析资料的时候是用numpy和pandas硬扛的写一堆for循环去读netCDF文件再手动记住哪个维度对应哪个轴换个数据源整套脚本就得推倒重来。后来真正系统地用上Xarray从读取、处理到可视化分析整条链路才算顺了起来。这一篇就把这套完整流程捋一遍先讲清楚Xarray的核心数据结构和工作方式再用一份真实的气象数据做读取、筛选、计算、重采样、插值最后用Cartopy把结果画成能直接放进报告和论文里的图。每一段都会给出能直接跑的Python代码并解释为什么这么写而不是只丢一个“标准答案”出来。适合刚被多维气象数据折磨过、想找一套可持续复用方案的读者也适合有pandas基础、打算进阶数据分析与可视化的人参考。1. 为什么是Xarray气象数据的天然结构与传统工具的痛点1.1 气象数据本质上是“带标签的多维数组”一说气象数据很多人第一反应是Excel表或者CSV但真实业务里接触最多的其实是netCDF、GRIB这类自描述格式。打开一份再分析资料里面往往是一个四维的数值场经度、纬度、时间、气压层或者高度层每个维度有自己的坐标值每个变量又带着单位、长名称等属性。比如一份温度数据结构可能是air(time8760, level17, lat73, lon144)。如果你是做气候统计的还要在这个基础上算月平均、季节平均、距平、区域平均。这种数据天生就不适合用“二维表格”的思维去处理因为一张表装不下这么多维度即便强行用MultiIndex撑起来代码也会变得又琐碎又难读。这也是我早期用pandas处理这类数据时最深的感受为了得到某一块区域的时间序列要先做切片、转置、reshape中间任何一步搞错维度顺序结果就全偏了而且很难排查。1.2 Xarray到底解决了什么问题Xarray的设计思路很直接在numpy多维数组的基础上给每个维度一个“名字”给每个轴一个“坐标标签”。这样所有操作都不再依赖“第0维是时间、第1维是纬度”这种脆弱的记忆而是直接通过名字和坐标值来操作。举个例子你不需要写data[2, 10:20, 30:40]这种只有自己看得懂的切片直接这样写就行data.sel(level850, latslice(20, 40), lonslice(100, 130))更关键的是Xarray有自动对齐机制。两个数据集的维度名字一致时做加减乘除会自动按坐标对齐不需要你手动对齐网格它还支持延迟计算配合Dask可以直接处理超出内存的大文件。这几点放在气象数据分析里几乎每一个都是刚需。还有一点容易被新手忽略Xarray的变量自带属性信息units、long_name处理完数据之后还能把这些信息原样写到结果文件里。后面你会发现这一点对流式复现和团队协作有多重要。2. 环境准备与数据源先把手上的工具磨利2.1 环境搭建优先用conda管理如果你还没有能用的Python环境我建议直接装一个Anaconda或者Miniconda。管理Python包这件事与其之后不断踩坑不如一开始就选一个省心的工具。气象领域常用的库很多依赖了编译好的底层库比如netCDF的C库、HDF5用conda装能省掉一大堆编译报错问题。我推荐单独建一个环境不要在base环境里堆一堆包conda create -n xarray_env python3.11 -y conda activate xarray_env然后安装本文会用到的库conda install -c conda-forge xarray netcdf4 dask python-cftime cartopy matplotlib -y这里几个依赖拆开解释一下netcdf4负责netCDF文件I/Opython-cftime用来处理带单位和参考日期的时间坐标cartopy是画地图投影和叠加海岸线的核心库dask则是给Xarray提供延迟计算和大数据分块能力。提示安装渠道强烈建议用conda-forge比默认的defaults渠道更新更快很多气象库都优先在conda-forge上发布。2.2 数据从哪来先拿真实数据练手如果你手头已经有实际的观测或模式输出nc文件那最好。如果没有可以直接用Xarray自带的教程数据也能把整套流程跑通import xarray as xr ds xr.tutorial.open_dataset(air_temperature) ds这份数据是一个经典的气温再分析场自带时间和经纬度坐标非常适合学习基本操作。等把流程跑通了再去下载真实的业务数据也不会手忙脚乱。真实业务中最常见的再分析资料比如ERA5、GFS、NCEP/NCAR这些都是公开可以获取的文件同样是netCDF或GRIB格式。拿到文件之后本文的读取和处理逻辑可以直接复用。唯一要提醒的是不同数据源的变量名、单位、坐标系可能有差异拿到新数据第一步先打印结构、确认attrs里的单位再开始动手分析。3. 核心数据结构与索引方式先看懂Dataset和DataArray3.1 Dataset是一个“容器”DataArray是里面的“数据块”刚开始接触Xarray的人最容易被这两个概念绕晕。我的理解方式是这样的Dataset可以看作一个文件夹里面可以放好几个变量DataArray就是文件夹里的单个文件是一个带坐标标签的多维数组。比如刚才打开的air_temperature它是一个Dataset里面只有一个变量air坐标是lat、lon、time。如果你想单独拿出来操作这个变量air ds[air] type(air) # xarray.core.dataarray.DataArrayDataArray有两个核心组成部分维度和坐标。维度是数组的轴名字比如(time, lat, lon)坐标是每个轴上具体的标签值比如lat坐标是[75.0, 72.5, ...]这类数值。理解这两者的区别很重要因为在做sel或者isel时一个按标签选一个按下标选用错了就会报错或者选中不是你要的数据。3.2 选择数据的四个常用方法sel、isel、loc、wheresel按坐标标签选取比如ds.sel(time2014-01-01)。isel按下标位置选取比如ds.isel(time0)表示取第一个时次。.loc在DataArray上也可以使用类似pandas的标签索引。where按条件筛选不符合条件的会变成NaN。这里有个很实用的点sel支持nearest选项当你需要找“离某一点最近的格点”时不必手动算索引位置# 找距离北京最近格点的温度序列 ds.sel(lat39.9, lon116.4, methodnearest)注意用sel做切片时slice的起止值是坐标值不是下标。如果坐标是降序排列比如纬度从90到-90latslice(20, 40)依然会正确切出20到40度之间的区间Xarray会自动处理顺序但如果你在操作前手动改过坐标排序一定要确认当前坐标是升序还是降序否则切出来的区域可能和你预想的不一样。3.3 用一个小例子快速建立直觉为了让自己对DataArray操作熟悉起来我建议先创建一个极简数据自己玩一玩import numpy as np import xarray as xr da xr.DataArray( np.random.rand(3, 4), dims(lat, lon), coords{ lat: [30, 32, 34], lon: [100, 102, 104, 106] }, nametemp, attrs{units: K} ) da.sel(lat32)这个例子虽然数据是随机数但你可以直观看到lat和lon是坐标temp是变量名单位信息也存在attrs里。后面所有复杂气象数据的处理本质上都是在这个小例子上做扩展。4. 完整处理流程从读取文件到输出结果4.1 读取文件open_dataset与open_mfdataset读取单个文件是open_dataset读多个文件用open_mfdataset。后者在气象数据分析里使用频率极高因为一整套再分析数据往往被拆成几年甚至几十年的多个文件无法一次性放到内存里。# 单个文件 ds xr.open_dataset(air_temperature.nc) # 多个文件自动按维度拼接 ds_multi xr.open_mfdataset(data/air_2014_*.nc, concat_dimtime, combinenested)open_mfdataset背后用的是Dask的延迟读取机制所以即便文件突然很大这行命令也不会把数据全部加载进内存只有后面真正需要计算的时候才会触发实际读取。这一点我建议从一开始就养成习惯大批量数据一律用open_mfdataset并合理设置chunks参数。ds xr.open_mfdataset( data/air_*.nc, concat_dimtime, combinenested, chunks{time: 100, lat: -1, lon: -1} )这里的chunks含义是每次处理100个时次纬度和经度整块加载。设置chunk大小的核心目标是让Dask的任务粒度适合你的计算资源不是越大越好也不是越小越好。实践上我一般先看单个文件的大小让每个chunk控制在几百MB以内再根据运行时的内存监控做调整。4.2 数据清洗缺失值、单位换算、区域裁剪拿到一份nc文件第一件事是检查结构、坐标范围、变量单位。我经常遇到的问题是变量的温度单位是开尔文画图时发现温度数值好几百那一定是没做单位换算。用属性信息来写通用处理逻辑是最稳的units ds[air].attrs.get(units, ) if units K: ds[air_c] ds[air] - 273.15 elif units degC: ds[air_c] ds[air] elif units degF: ds[air_c] (ds[air] - 32.0) * 5.0 / 9.0 else: # 默认按开尔文处理并打印提醒 ds[air_c] ds[air] - 273.15 print(未知单位按开尔文处理:, units)这样写的好处是即使下周换了一份数据源这段逻辑也能直接复用。缺失值在气象数据里通常以NaN形式存在。Xarray里的很多统计函数默认会跳过NaN比如mean(skipnaTrue)但在做加减乘除的时候NaN会“传染”。所以建议在关键计算之前先确认数据中有没有缺失值# 统计每个时次的NaN数量 missing ds[air_c].isnull().sum(dim(lat, lon))如果某个区域数据一直缺失那就老老实实做插值或者直接mask掉不要让NaN悄悄影响后面的距平和趋势计算。区域裁剪也是数据清洗的常见步骤用sel配合slice就能实现region ds[air_c].sel(latslice(20, 50), lonslice(100, 135))注意如果原始数据Lon是0到360而你习惯看-180到180的坐标先转换一下ds ds.assign_coords(lon(((ds[lon] 180) % 360) - 180)).sortby(lon)这个技巧在处理全球网格数据时几乎每次都能用上。4.3 重采样与聚合日均转月均、季节平均、滚动平均气象数据分析里重采样是走不掉的环节。再分析数据经常是一小时一次或六小时一次但做气候统计分析时通常需要日均、月均甚至季节平均。日平均到月平均最简单一行代码monthly ds[air_c].resample(time1M).mean()这里1M是“月结束”的频别名Xarray会按自然月分组然后对组内所有时次求平均。如果你想要“月始”对齐可以用MS两者差异在于时间戳标在月末还是月初。季节平均可以结合groupbyseasonal ds[air_c].groupby(time.season).mean()season是Xarray根据时间坐标自动识别出的季节标签返回结果是四个季节的平均场很好用。不过要提醒一下这里的季节是按气象季节划分的DJF、MAM、JJA、SON如果你想按1-12月的自然月平均用groupby(time.month)。还有滚动平均常用来做时间序列平滑smooth ds[air_c].rolling(time7, centerTrue).mean()这里的centerTrue表示滑动窗口以当前时次为中心适合低温滤波如果你需要的是未来时刻的信息不能用中心窗口就设置centerFalse。4.4 科学计算距平、梯度与通量估算很多人在这一步容易卡住因为气象里的计算往往不是一个简单公式就能做完的涉及单位换算、坐标重投影、差分方向等。Xarray的好处是它保留了维度和坐标信息差分计算变得非常直观。先看距平计算这是气候分析的重点。逐日数据通常需要减去“多年逐日气候态”得到的就是距平场climatology ds[air_c].groupby(time.month).mean(time) anomaly ds[air_c].groupby(time.month) - climatology这段代码的思路是先把所有年份同一个月的值做平均得到“这个月的多年平均态”再用原始值减去这个气候态得到每个时次相对当月的偏差。这是做温度异常、极端事件分析最常见的操作之一。如果要算温度的水平梯度Xarray也有直接的差分方法# 需要先确保lat/lon是等距的否则要用实际距离做差分 dT_dx ds[air_c].differentiate(lon) dT_dy ds[air_c].differentiate(lat)这里要注意气象网格的经度间隔和纬度间隔对应的实际距离不同如果要做严格意义上的梯度需要乘上cos(lat)修正。但在很多简单的分析任务里直接对坐标差分已经足够看到空间分布特征。如果你要做散度、涡度这类涉及二维差分和矢量分解的计算可以用xarray的apply_ufunc方法调用NumPy或者Scipy里的自定义函数。这个方法稍微有点学习成本但一旦理解就能处理任何逐格点运算。最基本的使用方式是这样的def compute_vorticity(u, v, dx, dy): dv_dx np.gradient(v, axis-1) / dx du_dy np.gradient(u, axis-2) / dy return dv_dx - du_dy vort xr.apply_ufunc( compute_vorticity, u, v, dx, dy, input_core_dims[[lon], [lon], [], []], output_core_dims[[]], daskparallelized, )核心要点是input_core_dims告诉Xarray哪些维度是变量的“核心维度”需要在传给函数之前把其他维度逐块迭代。一开始上手会觉得有些别扭但它是将普通NumPy函数升级成支持维度感知、延迟计算的标准做法。4.5 网格插值与坐标转换处理不同来源的数据时最烦的就是网格不一致。比如一个数据是1°分辨率另一个是0.25°分辨率要叠加对比时就必须插值到同一网格上。Xarray的interp做得非常顺手new_lat np.linspace(20, 50, 61) new_lon np.linspace(100, 135, 71) ds_interp ds[air_c].interp(latnew_lat, lonnew_lon)interp默认使用线性插值对大多数气象变量够用了。如果你是做模式输出对比建议用interp之后先对比一下结果的空间分布别盲信插值出来的数值特别是在地形复杂区域。还有一个常见的场景是把数据从规则经纬网格转换到某个特定站的单点序列。这个用前面的sel(methodnearest)就够了但如果站点和最近格点距离较远最好用双线性插值提取# 提取单点时间序列 station_ts ds[air_c].interp(lat39.9, lon116.4)插值有一个隐坑如果目标点落在源数据范围之外返回的会是NaN。比如源数据只覆盖到北纬20到50度你插值到北纬60度结果就是NaN。排查插值结果异常时先检查目标网格的范围是否在源网格的范围之内。4.6 结果保存与转pandas处理完的气象场我通常直接保存成netCDF文件这样变量属性和坐标信息都能完整保留下次读取还能继续用ds_interp.to_netcdf(air_interp_2014.nc)如果你想把结果交给只会用Excel的同事也可以转成pandas的DataFrame导出CSVdf ds_interp.isel(time0).to_dataframe() df.to_csv(air_2014-01-01.csv)不过我强烈建议只要是多维气象数据优先用netCDF保存。CSV会损失维度结构和属性信息等你下次再想画图、再想按时间切片又得重新整理一遍。5. 可视化把分析结果变成看得懂的图5.1 画图预备投影、海岸线、地图要素数据算完之后可视化是让结果“说话”的关键。Xarray自带的.plot方法可以用来快速预览但要做一张能放进正式报告里的地图我还是推荐直接使用matplotlib加cartopy。先设置一个基础地图画布import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.coastlines(linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.OCEAN, colorlightblue, alpha0.3) ax.set_extent([100, 135, 20, 50], crsccrs.PlateCarree())这里PlateCarree是等经纬度投影做中小尺度区域分析最常用。如果你画全球场建议换成ccrs.Robinson()或者ccrs.EqualEarth()视觉效果会好很多。5.2 填色图与等值线一张能放进报告的气温图画单时次的温度场最常用的是填色图叠加等值线。填色图能让空间分布一眼看清等值线则提供精确的数值参考import matplotlib.ticker as mticker data region.isel(time10) # 取第10个时次 fig plt.figure(figsize(12, 5)) ax plt.axes(projectionccrs.PlateCarree()) cf ax.contourf( data[lon], data[lat], data, levels20, cmapRdBu_r, transformccrs.PlateCarree() ) cs ax.contour( data[lon], data[lat], data, levels10, colorsk, linewidths0.8, transformccrs.PlateCarree() ) ax.clabel(cs, cs.levels[::2], fontsize8, fmt%.0f) ax.coastlines(linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5) ax.set_extent([100, 135, 20, 50], crsccrs.PlateCarree()) cbar plt.colorbar(cf, axax, orientationhorizontal, shrink0.7, pad0.05) cbar.set_label(Temperature ($^\circ$C)) ax.set_title(Air Temperature) plt.show()这段代码里最值得注意的地方是transformccrs.PlateCarree()。因为数据坐标是经纬度而画布坐标是投影坐标cartopy需要知道数据的坐标参照系统才能正确换算。这个参数在叠加多图层时尤其容易漏掉漏掉的后果是图形位置偏移肉眼还不容易发现。等值线标注那里我用cs.levels[::2]做抽样避免标注太密。你还可以用plt.xticks(rotation45)或者手动设置刻度间隔解决横坐标标签过密的问题。5.3 风场箭头图与风羽风场可视化通常有两种方式箭头图quiver和风羽图barbs。区域尺度的风场用箭头更直观展示风向和相对风速大小天气图风格则常用风羽。fig plt.figure(figsize(12, 5)) ax plt.axes(projectionccrs.PlateCarree()) # 假设u10、v10是10米风场先做抽样避免箭头太密 uq u10.isel(time0).sel(latslice(20, 50), lonslice(100, 135)) vq v10.isel(time0).sel(latslice(20, 50), lonslice(100, 135)) ax.quiver( uq[lon][::4], uq[lat][::4], uq[::4, ::4], vq[::4, ::4], transformccrs.PlateCarree(), scale300, width0.002 )scale参数控制箭头长度不同区域风速默认表现差异很大需要手动调。新手最容易犯的错是把箭头画得太密或太长视觉上一团乱麻。我的习惯是先抽样[::4]再打开图看一次根据最大风速值调整scale。风羽图适合展示风速等级和风向适合做更专业的天气分析ax.barbs( uq[lon][::3], uq[lat][::3], uq[::3, ::3], vq[::3, ::3], transformccrs.PlateCarree(), length5, barb_incrementsdict(half5, full10, flag50) )很多气象观测站的台风、大风报告图用的就是这种风羽风格。5.4 时间序列、剖面与多子图组合区域平均后的时间序列图是展示一个区域温度演变最直观的方法regional_mean region.mean(dim(lat, lon)) regional_mean.plot(figsize(10, 3)) plt.ylabel(Temperature ($^\circ$C))这里可以看到Xarray的.plot针对DataArray的默认行为一维数据自动画折线图。如果横坐标时间点太多导致标签拥挤可以用plt.gca().xaxis.set_major_locator(plt.MaxNLocator(6))强制只显示少量刻度避免一坨标签叠在一起。剖面图也是气象分析的高频需求。假设你的数据有level维可以用等值线填色画纬向垂直剖面fig, ax plt.subplots(figsize(10, 4)) prof temp.mean(lon) # 经向平均后剩下 (level, lat) cf ax.contourf(prof[lat], prof[level], prof, levels20, cmapviridis) ax.invert_yaxis() # 气压层从大到小y轴翻转后高层在上 ax.set_yscale(log) # 气压层常用对数坐标关于气压层坐标这里有个细节气压层的数值从地面到高空是递减的大气科学里通常把y轴设置成“上小下大”因此要invert_yaxis()。如果还用线性坐标近地面变化大的区域会被压缩得看不清所以常用对数或自定义坐标。多子图组合比如把填色图和时间序列放一起fig, (ax1, ax2) plt.subplots( 1, 2, figsize(15, 5), subplot_kw{projection: ccrs.PlateCarree()} )第一个子图画空间分布第二个子图画该区域的均值时间序列。这种组合图在写总结报告时非常实用一页就能讲清“空间上哪里偏暖、时间上怎么演变”。6. 常见问题与排查技巧实录6.1 踩过的坑速查表在整理这套流程的这段日子里我确实积累了一些经典问题下面这个表格基本覆盖了新手最常见的情况。报错现象可能原因排查与解决KeyError: time或时间索引报错缺少cftime依赖或时间坐标不是标准datetime安装python-cftime用ds.indexes[time]检查类型resample报“time must be monotonic”时间坐标有重复或乱序先ds ds.sortby(time)再用drop_duplicates(time)插值结果大量NaN目标网格超出源数据坐标范围打印源坐标范围[min, max]确认目标点落在范围内画图位置整体偏移少了transformccrs.PlateCarree()所有用经纬度坐标绘制的图层都要声明transform内存占用异常高数据被全部load进内存改用open_dataset(..., chunks...)或open_mfdataset计算时保持延迟groupby(time.month)算出来全是0或NaN时间坐标不是datetimelike类型先转pd.to_datetime并赋给时间坐标区域切片结果反了坐标是降序slice理解错了先用ds.lat和ds.lon打印范围再用sortby保证升序to_dataframe()后行数爆炸Dataset里有多个不同shape的变量先选择你要的变量再转避免坐标笛卡尔积其中“时间坐标问题”是真正的重灾区。很多nc文件的时间单位是hours since 1900-01-01这种CF标准格式Xarray需要cftime才能正确解码。如果你不想每次都被这个坑绊倒建议在环境创建时直接把python-cftime装上。6.2 几个能提升幸福感的小习惯有些习惯看起来不起眼但在长期使用Xarray之后它们真的能让你少走很多弯路。第一写处理函数时尽量接受DataArray或Dataset作为输入返回同类型对象。这样一来所有函数都能用.pipe串起来整个流程变成一条清晰的流水线result ( ds .pipe(decode_units) .pipe(resample_monthly) .pipe(calculate_anomaly) )第二处理完结果后把关键信息写进attrs。Xarray允许你在变量上直接挂属性这些属性会在保存netCDF时一并写入。以后翻结果文件的时候不需要再回到原始脚本里去猜“这个变量当时是怎么算的”。result.attrs[description] 每月平均的2米气温距平 result.attrs[history] 由raw_air_temperature_2014.nc计算所得 result.attrs[units] degC第三别急着用大数组调试。我在整合大文件前习惯先对coords做一次抽样比如只取两三个时次、小块区域跑通再放大到全量数据。这个小习惯帮我节省了大量等待时间尤其是第一次跑不熟悉的计算逻辑时。第四多用.plot()快速检查中间结果。Xarray的DataArray自带绘图方法一维数据自动折线二维自动填色四维数据需要先isel或sel到二维。做复杂处理之前先画图看分布是否合理比闷头算一个方差、相关系数更有效。第五关于apply_ufunc如果确实需要写自定义的逐格点运算建议先把小块的NumPy函数测试好再包给Xarray。不要在一个巨大数据集上直接调试否则一次报错就是几分钟的时间成本。最后再分享一点我自己的体会从开始用Xarray到现在我觉得最值得的投入其实不是学各种API而是理解“维度感知”这一种思维模式。过去用numpy处理数据我脑子里永远绷着一根弦某个轴顺序不能记反。现在用Xarray我只需要关心变量名和坐标名剩下的一切都由库来管理。我的工作流里现在几乎每条数据处理路径都用Xarray打底哪怕只是做一个小范围的统计我也会先把nc文件读成Dataset再说。这套流程我建议你不要一口气全部跑完而是按章节拆开跑。先读懂数据结构和选择方法再动手做重采样和计算最后画图。每一个环节都跑通之后把代码整理成自己的函数库后面处理任何再分析资料都会快很多。如果你刚开始接触多维气象数据第一个小目标就定成把一份nc文件读进来截图看一下数据结构再试着画一张和本文类似的图。跑通这一步后面就是水到渠成的事。