pangeo xarray ENSO练习:SST异常计算与ONI绘图问题求解
已跑通的基础流程
- 克隆配套数据仓库的命令:
git clone https://github.com/pangeo-data/tutorial-data.git
- 基础数据加载与预处理代码:
import numpy as np import xarray as xr import pandas as pd # 加载全部SST数据文件 ds_all = xr.open_mfdataset('./tutorial-data/sst/*nc', combine='by_coords') # 截取Niño3.4研究区域 sst_enso = ds_all.sel(lat=slice(-5,5), lon=slice(-170+360,-120+360)) # 截取1982年之后的时间范围 sst_enso = sst_enso.sel(time=sst_enso.time.dt.year>=1982) # 计算区域平均SST时间序列 avg_enso = sst_enso.sst.groupby('time').mean(dim=['lat','lon']) avg_enso.plot()
这部分流程输出的Niño3.4区域平均SST图与教程示例效果一致,可正常推进后续计算。
核心操作卡点
卡点1:SST异常计算的维度匹配疑问
尝试按月分组计算SST气候态并计算异常,测试代码为:
enso_clim = sst_enso.sst.groupby('time.month').mean('time') sst_anom = sst_enso - enso_clim enenso_anom = sst_anom.groupby('time').mean(dim=['lat','lon']) oni = enso_anom.rolling(time = 3).mean()
因groupby操作引入了month维度索引,无法确认原始SST场与气候态场的逐元素相减逻辑是否正确,担心异常计算存在维度匹配偏差。
卡点2:冗余month维度导致时间序列异常
完成ONI序列计算后,尝试绘制正红负蓝的分色柱状图,测试代码为:
oni.sst.plot.bar(color=(oni.sst < 0).map({True: 'b', False: 'r'}))
实际运行出现四类异常:
- 普通线图输出为年-月交叉的二维网格图,而非预期的单时间序列
- 用
reset_index('month', drop=True)重置索引后,month维度仍未移除,绘图效果无变化 - 直接调用
drop_dims('month')会丢失原有SST数据 - 将oni转换为pandas DataFrame时得到5040行冗余结果,即12个月份×420个年-月时间点。根据xarray的坐标笛卡尔积构建逻辑,即使转换前重置oni索引,冗余结果也不会消失,无法得到无冗余的单时间步序列完成分色柱状图绘制。
对应解决方法
解决SST异常匹配问题
xarray的groupby运算自带自动对齐匹配机制,只要按time.month分组得到气候态,和原始场做减法时,系统会自动把每个时间点的SST值和对应月份的气候态值做匹配相减,不需要手动做维度对齐操作。
之前测试代码里的sst_anom = sst_enso - enso_clim写法,是用包含经纬度、时间维度的完整数据集,和单独的(month, lat, lon)三维气候态数组做广播运算,系统会把气候态带的month维度单独识别为一个新的独立匹配维度,这就是后续出现month维度残留、生成二维年-月网格图的根本原因。
正确的异常计算流程应把减法嵌套在同个groupby逻辑里完成,从根源上避免引入额外维度:
# 计算多年逐月气候态 enso_clim = sst_enso.sst.groupby('time.month').mean('time') # 同groupby流程内做减法,自动匹配对应月份的气候态 sst_anom = sst_enso.sst.groupby('time.month') - enso_clim # 空间平均得到逐时间点的区域异常序列 enso_anom = sst_anom.mean(dim=['lat', 'lon']) # 3个月滑动平均得到ONI序列 oni = enso_anom.rolling(time=3).mean()
解决冗余month维度问题
按照上述规范流程运算后,enso_anom和oni本身就只保留单时间维度,不会额外引入month维度,自然不会出现笛卡尔积的冗余结果。很多人遇到的转换后5040行冗余问题,本质就是之前广播运算时引入了多余的month维度——xarray转pandas时会对所有坐标做笛卡尔积展开,12个month值乘上420个time值,就会得到12倍的冗余行数。
如果运算后确实残留了不需要的month坐标,不要用容易丢数据的drop_dims方法,用reset_coords即可安全清理:
# 安全移除残留的month坐标,不影响实际数据值 oni = oni.reset_coords('month', drop=True)
分色柱状图的正确绘制逻辑需要传入和时间长度完全一致的颜色列表,不要直接对可能为二维的布尔序列做map运算,正确写法为:
import matplotlib.pyplot as plt # 逐值分配颜色:负值蓝色,正值红色 colors = np.where(oni.values < 0, 'b', 'r') # 绘制ONI柱状图 oni.plot.bar(color=colors, width=0.8) # 添加±0.5℃的阈值参考线,标记厄尔尼诺/拉尼娜的判定标准 plt.axhline(y=0.5, c='r', ls='--', lw=0.8) plt.axhline(y=-0.5, c='b', ls='--', lw=0.8)
这种方法得到的结果不会出现冗余的笛卡尔积匹配,每个时间点只会对应一个颜色值,输出的单序列柱状图符合预期。
内容的提问来源于stack exchange,提问作者Dominik

