Python气象数据处理进阶之Xarray(7)：读写文件

Posted 2023-05-02

参考技术A 前几文主要讲的是如何处理Xarray中的DataArray和DataSet，现在分享一下如何从nc文件或其他文件中读取数据，以及如何将处理好的数据输出成一个nc文件。

首先还是要再强调DataArray和DataSet的区别，DataArray是一个带标签结构的数组，DataSet是一个数据集，这意味着，从一个nc文件中读取到的全部信息构成了一个DataSet，而nc文件中的某一个变量是一个DataArray。

反之，我们要将一个数据写成nc文件，那么就是要创建一个DataSet。

这个数据结构有点像站点数据，对xy维设定了两层，分别是经纬度，还有一维时间维(whatever，反正是随便创建一个DataSet)。

就可以输出成nc文件了。
当然还可以更懒一点，

直接将abc这个DataArray转成DataSet，DataArray的标签和纬度信息会自动转换。
之后使用to_netcdf即可。

读取的语句也十分简单。

函数只需要基本的路径及文件名，无需像NCL一样声明状态'r'。

Xarray读取多文件也提供了相应函数(我目前没有使用过，我通常都是使用CDO提前处理，大家可以自行尝试)。

根据官方的介绍，Xarray也支持grib文件的读取。

前提是需要一个解码库"eccodes"

或者利用Xarray借助PYNIO去读。

官方文档中还有一部分是关于画图的，然而画图部分个人认为使用matplotlib+cartopy的组合更加灵活，因此Xarray系列到这里应该就完结了。

下一步的计划是按照魏凤英老师的统计方法一书，试着将常用的气象统计方法利用python去实现，但是水平实在有限。

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

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

【英文标题】：How to join data from multiple netCDF files with xarray in Python?

【发布时间】：2019-08-22 18:45:55

【问题描述】：

我正在尝试在 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）。 它创建了一个新维度，但也将最左边的维度相加：两个数据集的“时间”维度。 我想知道这是来自“xarray.open_mfdateset”的错误还是预期的行为，因为两个数据集的“时间”维度都是无限的。

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

谢谢。

马特乌斯

【问题讨论】：

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

两个文件的时间维度长度相同，124个时间索引。但是，它们不涵盖相同的时间段。

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

如果您考虑具有不同时间段（但长度相等）的不同案例或事件，您可能有一个包含不同事件数据的大型数组，每个事件都有自己的持续时间。因此，结果数组将有一个额外的维度：事件。

@Mateus 您可以尝试使用每个文件都相等的新维度（仅 0,1，...，123）来重命名/复制/删除时间维度。这可以通过 open_mfdataset 的 preprocess 关键字参数来完成。我正在做类似的事情，但手头上没有代码。

【参考方案1】：

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

为了简化它，暂时忘记 lat-lon 维度，假设您有两个文件，它们只是 2 个时间片上的数据。第一个文件的时间步长为 1,2，第二个文件的时间步长为 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 文件，时间维度必须跨越两个文件中存在的值范围。您可以合并两个文件并获取（案例：2，时间：124，纬度：241，经度：480）的唯一方法是如果两个文件具有相同的时间、纬度和经度值，即指向完全相同的区域时间-纬度空间。

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

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

【参考方案2】：

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

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

我正在尝试做与 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

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

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

再次感谢大家的帮助！

马特乌斯

【讨论】：

我明白了，为什么不在循环中分别读取每个文件并将数据放入 4D numpy 数组中，而不是使用 xarray concat 函数？

我只是想避免循环，使用这个 xarray 功能。但是，我想我做不到。谢谢。

【参考方案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)

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

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

免责声明：我在最后一个 concat 的循环中做了类似的事情（效果很好），但没有测试preprocess-方法。

【讨论】：

谢谢@kmuehlbauer！您的代码运行良好，包括将日期保存在 stime 变量中。

向所有在 PyData 保护伞下做出出色工作的人致敬。这些包（例如 pandas、xarray 等）确实改善了我的工作流程。