减去两个 xarray，同时保留所有维度

Question

这可能是最基本的问题，但我就是找不到解决方案。

我有两个包含风数据的不同 xarray。两个 xarray 都有维度（时间：60，普列夫：19，纬度：90）。我现在需要在所有维度上获取两个 xarray 之间的差异，以找到两个场景之间的异常。

我认为 xarray.DataArray.diff 函数仅用于计算沿一个 xarray 的轴的差异（而不是计算两个 xarray 之间的差异）。

所以，我尝试简单地使用

diff = wind1_xarray - wind2_xarray

也

diff = (wind1_xarray - wind2_xarray).compute() 

然而，这两种方法都给了我一个具有尺寸的xarray（时间：60，plev：0，lat：90）。为什么在计算差异时会丢失压力水平？

如何计算两个 xarray 在所有维度上的差异而不丢失一维？

感谢大家

Answer 1

简单的回答是，你做得对，但你的维度没有对齐。xarrayIS 设计用于减去整个数组，但坐标标签必须精确对齐。您的 plev 坐标元素之间可能存在分歧，您可以通过以下方式进行检查xr.align：

xr.align(wind1_array, wind2_array, join='exact')

有关更多信息，请参阅有关计算的 xarray 文档：自动对齐。

详细示例

xarray 和 numpy 之间的最大区别（假设您熟悉使用 numpy 进行数学计算）是 xarray 依赖于每个维度上的坐标标签来在任何广播操作之前对齐数组，而不仅仅是形状。

作为示例，让我们考虑两个非常简单的数组 - 一个从 0 计数到 19，另一个是一组 1，两者都重新整形为(4, 5)。在 numpy 中将它们相减很简单，因为它们的形状相同：

In [15]: arr1 = np.arange(20).reshape((4, 5))

In [16]: arr2 = np.ones(shape=(4, 5))

In [17]: arr1 - arr2
Out[17]:
array([[-1.,  0.,  1.,  2.,  3.],
       [ 4.,  5.,  6.,  7.,  8.],
       [ 9., 10., 11., 12., 13.],
       [14., 15., 16., 17., 18.]])

xarray 等效项也很简单，但我们必须引入维度名称和坐标。假设您的气压水平以 10 hPa 的增量向 STP 递减，纬度也以 10 的增量从 20 到 60：

In [18]: pressures = np.array([71.325, 81.325, 91.325, 101.325])

In [19]: lats = np.array([20, 30, 40, 50, 60])

In [20]: da1 = xr.DataArray(arr1, dims=['plev', 'lat'], coords=[pressures, lats])

In [21]: da2 = xr.DataArray(arr2, dims=['plev', 'lat'], coords=[pressures, lats])

In [22]: da2
Out[22]:
<xarray.DataArray (plev: 4, lat: 5)>
array([[1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.]])
Coordinates:
  * plev     (plev) float64 71.33 81.33 91.33 101.3
  * lat      (lat) int64 20 30 40 50 60

In [23]: da1
Out[23]:
<xarray.DataArray (plev: 4, lat: 5)>
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])
Coordinates:
  * plev     (plev) float64 71.33 81.33 91.33 101.3
  * lat      (lat) int64 20 30 40 50 60

这些数组是对齐的，因此减去它们很简单：

In [24]: da1 - da2
Out[24]:
<xarray.DataArray (plev: 4, lat: 5)>
array([[-1.,  0.,  1.,  2.,  3.],
       [ 4.,  5.,  6.,  7.,  8.],
       [ 9., 10., 11., 12., 13.],
       [14., 15., 16., 17., 18.]])
Coordinates:
  * plev     (plev) float64 71.33 81.33 91.33 101.3
  * lat      (lat) int64 20 30 40 50 60

但由于 xarray 依赖于这些坐标的精确对齐，因此依赖浮点坐标可能会很棘手。如果我们在压力水平维度上引入哪怕一个小误差，数组就会不对齐，我们会看到与您类似的结果：

In [25]: da2 = xr.DataArray(arr2, dims=['plev', 'lat'], coords=[pressures + 1e-8, lats])

In [26]: da1 - da2
Out[26]:
<xarray.DataArray (plev: 0, lat: 5)>
array([], shape=(0, 5), dtype=float64)
Coordinates:
  * plev     (plev) float64
  * lat      (lat) int64 20 30 40 50 60

这种类型的未对齐可能因各种原因而发生，包括通过存储来回传输数据，其中编码的更改可能会导致微小的数字错误，这些错误会显示为未对齐的数据。

xr.align您可以检查这是否是使用参数的问题join='exact'：

In [27]: xr.align(da1, da2, join='exact')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-29-612460e52308> in <module>
----> 1 xr.align(da1, da2, join='exact')

~/miniconda3/envs/myenv/lib/python3.9/site-packages/xarray/core/alignment.py in align(join, copy, indexes, exclude, fill_value, *objects)
    320             ):
    321                 if join == "exact":
--> 322                     raise ValueError(f"indexes along dimension {dim!r} are not equal")
    323                 joiner = _get_joiner(join, type(matching_indexes[0]))
    324                 index = joiner(matching_indexes)

ValueError: indexes along dimension 'plev' are not equal

要解决此问题，您可以尝试将坐标舍入为已知的坐标容差：

In [32]: da2['plev'] = np.round(da2['plev'], 3)

In [33]: da1 - da2
Out[33]:
<xarray.DataArray (plev: 4, lat: 5)>
array([[-1.,  0.,  1.,  2.,  3.],
       [ 4.,  5.,  6.,  7.,  8.],
       [ 9., 10., 11., 12., 13.],
       [14., 15., 16., 17., 18.]])
Coordinates:
  * plev     (plev) float64 71.33 81.33 91.33 101.3
  * lat      (lat) int64 20 30 40 50 60

或者，您可以设置位置/整数坐标，将实际压力水平作为非索引坐标：

In [42]: da1
Out[42]:
<xarray.DataArray (plev_ind: 4, lat: 5)>
array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19]])
Coordinates:
    plev      (plev_ind) float64 71.33 81.33 91.33 101.3
  * lat       (lat) int64 20 30 40 50 60
  * plev_ind  (plev_ind) int64 71325 81325 91325 101325

In [43]: da2
Out[43]:
<xarray.DataArray (plev_ind: 4, lat: 5)>
array([[1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1.]])
Coordinates:
    plev      (plev_ind) float64 71.33 81.33 91.33 101.3
  * lat       (lat) int64 20 30 40 50 60
  * plev_ind  (plev_ind) int64 71325 81325 91325 101325

In [44]: da1 - da2
Out[44]:
<xarray.DataArray (plev_ind: 4, lat: 5)>
array([[-1.,  0.,  1.,  2.,  3.],
       [ 4.,  5.,  6.,  7.,  8.],
       [ 9., 10., 11., 12., 13.],
       [14., 15., 16., 17., 18.]])
Coordinates:
  * lat       (lat) int64 20 30 40 50 60
  * plev_ind  (plev_ind) int64 71325 81325 91325 101325