使用 xarray 滚动的有限差分

Question

我的目标是计算多维数据集沿给定维度的移动窗口的导数，其中数据集存储为 Xarray DataArray 或 DataSet。

在最简单的情况下，给定一个二维数组，我想计算一维中多个条目的移动差异，例如：

data = np.kron(np.linspace(0,1,10), np.linspace(1,4,6) ).reshape(10,6)
T=3
reducedArray = np.zeros_like(data)
for i in range(data.shape[1]):
    if i < T:
        reducedArray[:,i] = data[:,i] - data[:,0]
    else:
        reducedArray[:,i] = data[:,i] - data[:,i-T]

if i <T 条件确保输入和输出包含正确的值（即，没有 nans）并且具有相同的形状。

Xarray 的diff 旨在使用最近邻对给定的导数阶进行有限差分逼近，因此这里不适合，因此问题是：

是否可以仅使用 Xarray 函数来执行此操作？

rolling weighted average 示例看起来有点相似，但由于使用了 NumPy 例程，仍然过于不同。我一直在想，应该按照以下思路进行操作：

xr2DDataArray = xr.DataArray(
    data,
    dims=('x','y'),
    coords={'x':np.linspace(0,1,10), 'y':np.linspace(1,4,6)}
)
r = xr2DDataArray.rolling(x=T,min_periods=2)
r.reduce( redFn )

不过，我在这里为redFn 的定义苦苦挣扎。

警告要应用该操作的实际数据集的大小约为 10GiB，因此我们将非常感谢您提供不会破坏内存需求的解决方案！

---

更新/解决方案

使用 Xarray rolling

睡在上面并稍微摆弄一下上面链接的帖子实际上包含一个解决方案。为了获得有限差分，我们只需将权重定义为在末端为 $\pm 1$，否则为 $0$：

def fdMovingWindow(data, **kwargs):
    T = kwargs['T'];
    del kwargs['T'];
    weights = np.zeros(T)
    weights[0] = -1
    weights[-1] = 1
    axis = kwargs['axis']
    if data.shape[axis] == T:
        return np.sum(data * weights, **kwargs)
    else:
        return 0

r.reduce(fdMovingWindow, T=4)

或者，使用construct 和点积：

weights = np.zeros(T)
weights[0] = -1
weights[-1] = 1
xrWeights = xr.DataArray(weights, dims=['window'])
xr2DDataArray.rolling(y=T,min_periods=1).construct('window').dot(xrWeights)

这带有一个大量警告：该过程实质上创建了一个表示移动窗口的列表数组。这对于一个普通的 2D / 3D 阵列来说很好，但是对于一个占用大约 10 GiB 内存的 4D 阵列来说，这将导致 OOM 死亡！

简单 - 内存高效

一种占用内存较少的方法是复制数组并以类似于 NumPy 的数组的方式工作：

xrDiffArray = xr2DDataArray.copy()
dy = xr2DDataArray.y.values[1]  - xr2DDataArray.y.values[0] #equidistant sampling
for src in xr2DDataArray:
    if src.y.values < xr2DDataArray.y.values[0] + T*dy:
        xrDiffArray.loc[dict(y = src.y.values)] = src.values - xr2DDataArray.values[0]
    else:
        xrDiffArray.loc[dict(y = src.y.values)] = src.values - xr2DDataArray.sel(y = src.y.values - dy*T).values

这将产生没有尺寸误差的预期结果，但它需要数据集的副本。

我希望利用 Xarray 来防止复制，而只是链接操作，然后在 如果和何时实际请求值时进行评估。

仍然欢迎提出如何完成此任务的建议！

Answer 1

我从来没有使用过xarray，所以也许我弄错了，但我认为你可以得到你想要的结果，避免使用循环和条件。这比 numpy 数组的示例至少快两倍：

data = np.kron(np.linspace(0,1,10), np.linspace(1,4,6)).reshape(10,6)
reducedArray = np.empty_like(data)
reducedArray[:, T:] = data[:, T:] - data[:, :-T]
reducedArray[:, :T] = data[:, :T] - data[:, 0, np.newaxis]

我想使用DataArrays 时改进会更高。 它不使用xarray 函数，但也不依赖numpy 函数。我相信将其翻译为xarray 会很简单，我知道如果没有coords，它会起作用，但是一旦包含它们，由于coords 不匹配（@ 987654330@ 的@ 987654331@ 和 data[:, :-T] 不同）。可悲的是，我现在不能做得更好。