如何使用 Python 中的 xarray 连接来自多个 netCDF 文件的数据？

Question

我正在尝试在 Python 中使用 xarray 打开多个 netCDF 文件。这些文件具有相同形状的数据，我想加入它们，创建一个新的维度。

我尝试对 xarray.open_mfdataset() 使用 concat_dim 参数，但它没有按预期工作。下面给出一个例子，打开两个温度数据文件124次，241个纬度和480个经度：

DS = xr.open_mfdataset( 'eraINTERIM_t2m_*.nc', concat_dim='cases' )
da_t2m = DS.t2m

print( da_t2m )

使用此代码，我希望结果数据数组的形状类似于（案例：2，时间：124，纬度：241，经度：480）。但是，它的形状是（案例：2，时间：248，纬度：241，经度：480）。它创建了一个新维度，但也对最左边的维度求和：两个数据集的 'time' 维度。我想知道这是 'xarray.open_mfdateset' 的错误还是预期的行为，因为两个数据集的 'time' 维度都是无限的。

有没有办法直接使用 xarray 连接这些文件中的数据并获得上述预期 return？

谢谢。

马特乌斯

Answer 1

如果时间不同，结果是有意义的。

为了简化它，暂时忘掉经纬度维度，假设您有两个文件，它们只是 2 个时间片的数据。第一个文件的时间为 timesteps 1,2，第二个文件的 timesteps 为 3 和 4。您不能创建时间维度仅跨越 2 个时间片的组合数据集；时间维度变量必须有时间 1,2,3,4。所以如果你说你想要一个新的维度 "cases"，然后数据被组合成一个二维数组，看起来像这样：

times: 1,2,3,4

cases: 1,2

data: 
               time
          1    2    3    4
cases 1:  x1   x2 
      2:            x3   x4

考虑等效的 netcdf 文件，时间维度必须跨越两个文件中存在的值范围。合并两个文件并获得 (cases: 2, time: 124, latitude: 241, longitude: 480) 的唯一方法是，如果两个文件具有相同的时间、纬度和经度值，即指向完全相同的区域时间经纬度 space.

ps：这个问题有点离题，但如果您刚刚开始新的分析，为什么不切换到新一代、更高分辨率的 ERA-5 再分析，现在可以返回也到 1979 年（最终会向后扩展），您可以使用此处的 python api 脚本将其直接下载到您的桌面：

https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset

Answer 2

谢谢@AdrianTompkins 和@jhamman。在您的评论之后，我意识到由于不同的时间段，我真的无法使用 xarray 得到我想要的东西。

我创建这样的数组的主要目的是在一个 N 维数组中获取具有相同持续时间的不同事件的所有数据。因此，我可以很容易地获得，例如，每个时间（小时、天等）所有事件的复合字段。

我正在尝试做与 NCL 相同的事情。请参阅下面的 NCL 代码，它对相同的数据（对我而言）按预期工作：

f = addfiles( (/"eraINTERIM_t2m_201812.nc", "eraINTERIM_t2m_201901.nc"/), "r" )
ListSetType( f, "join" )
temp = f[:]->t2m
printVarSummary( temp )

最终结果是一个4维数组，新的自动命名为ncl_join.

但是，NCL 不考虑时间轴，连接数组并将第一个文件的坐标提供给生成的时间轴。所以，时间轴就没用了。

但是，正如@AdrianTompkins 所说，时间段不同，xarray 无法像这样连接数据。因此，要创建这样的数组，在 Python 中使用 xarray，我认为唯一的方法是从数组中删除时间坐标。因此，时间维度将只有整数索引。

xarray 给出的数组就像@AdrianThompkins 在他的小例子中所说的那样工作。由于它保留所有合并数据的时间坐标，与 NCL 相比，我认为 xarray 解决方案是正确的解决方案。但是，现在我认为复合材料的计算（得到上面给出的相同示例）不会像 NCL 看起来那么容易。

在一个小测试中，我用

打印了来自合并数组的两个值和 xarray

print( da_t2m[ 0, 0, 0, 0 ].values )
print( da_t2m[ 1, 0, 0, 0 ].values )

结果是什么

252.11412
nan

第二种情况，第一次没有数据，符合预期。

UPDATE：所有答案都帮助我更好地理解这个问题，所以我不得不在这里添加一个更新，也感谢@kmuehlbauer 的回答，表明他的代码给出了预期数组。

再次感谢大家的帮助！

马特乌斯

Answer 3

从我的评论延伸我会试试这个：

def preproc(ds):
    ds = ds.assign({'stime': (['time'], ds.time)}).drop('time').rename({'time': 'ntime'})
    # we might need to tweak this a bit further, depending on the actual data layout
    return ds

DS = xr.open_mfdataset( 'eraINTERIM_t2m_*.nc', concat_dim='cases', preprocess=preproc)

这里的好处是，在重命名原始维度 (time -> ntime) 的同时，您在 stime 中保留了原始时间坐标。

如果一切正常，您应该得到结果尺寸为 (cases, ntime, latitude, longitude).

免责声明：我在带有最终连接的循环中做了类似的事情（效果很好），但没有测试 preprocess 方法。

如何使用 Python 中的 xarray 连接来自多个 netCDF 文件的数据？

How to join data from multiple netCDF files with xarray in Python?

python

concatenation

netcdf

python-xarray