如何使用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 )

使用这段代码，我希望结果数据数组的形状为 (cases: 2, time: 124, latitude: 241, longitude: 480)。然而，实际上它的形状是 (cases: 2, time: 248, latitude: 241, longitude: 480)。这个代码创建了一个新的维度，但也对左侧维度（即两个数据集的 'time' 维度）进行了求和。

我在想，这是由于 'xarray.open_mfdateset' 的错误还是因为两个数据集的 'time' 维度都是无限的而导致的预期行为？

有没有一种方法可以直接使用 xarray 将这些文件中的数据连接起来，并获得所期望的结果？

谢谢。

Mateus

- Mateus da Silva Teixeira

---

时间维度是否相同（即两个文件是否包含相同时间范围的数据）？ - user7813790

时间维度在两个文件中具有相同的长度，共有124个时间索引。然而，它们并未涵盖相同的时间段。 - Mateus da Silva Teixeira

Xarray 正在尝试沿着现有的时间维度对齐您的数据。在 xarray 连接各个数据集之后，您希望您的时间坐标是什么？ - jhamman

如果您考虑到具有不同时间段（但长度相等）的不同情况或事件，您可能会拥有一个包含不同事件数据的大型数组，每个事件都有自己的持续时间。因此，结果数组将具有一个附加维度：事件。 - Mateus da Silva Teixeira

@Mateus 你可以尝试通过open_mfdataset的preprocess关键字参数，将时间维度重命名/复制/删除为一个新的维度，该维度在每个文件中都相等（只需为0,1,……，123）。我正在做类似的事情，但此时手头没有代码。 - kmuehlbauer

Answer 1

扩展我的评论，我会尝试这样做：

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方法。

- kmuehlbauer

---

谢谢@kmuehlbauer！你的代码非常好用，包括将日期保存在stime变量中。 - Mateus da Silva Teixeira

向所有在 PyData 傘下做出卓越工作的人致敬。这些软件包（例如 pandas，xarray 等）真正改进了我的工作流程。 - kmuehlbauer

Answer 2

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

为了简化，暂时忘记经纬度维度，并想象你有两个文件，它们只是在2个时间片段上的数据。第一个文件在时间步长1、2处有数据，第二个文件在时间步长3和4处有数据。您无法创建一个时间维度仅跨越2个时间片段的组合数据集；时间维度变量必须具有1、2、3、4这些时间。因此，如果您说要一个新的维度“cases”，那么数据将被组合为2D数组，并且看起来像这样：

times: 1,2,3,4

cases: 1,2

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

想象一下相应的netcdf文件，时间维度必须跨越两个文件中存在的值范围。只有当两个文件具有相同的时间、纬度和经度值时，才能合并两个文件并获得（案例：2，时间：124，纬度：241，经度：480）。也就是说，它们指向时间-纬度-经度空间中完全相同的区域。

附注：这与问题有些偏题，但如果您刚开始进行新的分析，为什么不转而使用新一代更高分辨率的ERA-5再分析数据呢？它现在也可用于1979年以来的数据（最终将会扩展到更久远），您可以从这里使用python api脚本直接下载到您的桌面：

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

- ClimateUnboxed

Answer 3

感谢 @AdrianTompkins 和 @jhamman 的评论。经过你们的指导，我意识到由于不同时间段的数据，我无法使用 xarray 达到我的目的。

我的主要目的是创建这样一个数组，在单个 N-D 数组中获取不同事件的所有数据，并且具有相同的时间跨度。因此，例如，我可以轻松地获取每个时间点（小时，天等）的所有事件的复合场。

我正在尝试与 NCL 相同的方法。以下是一段在 NCL 中可以正常工作（对我而言）的代码：

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

最终结果是一个四维数组，新的数组自动命名为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

对于第二种情况，预期的第一次没有数据。

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

再次感谢大家的帮助！

Mateus

- Mateus da Silva Teixeira

---

我明白了，为什么不在循环中分别读取每个文件，并将数据放入一个4D numpy数组中呢？这样做比使用xarray concat函数更好。 - ClimateUnboxed

我只是想避免使用循环，利用xarray的能力。但是，我认为我可能做不到。谢谢。 - Mateus da Silva Teixeira