python - 跨时间暗淡的每个网格单元的线性回归

Question

我是 xarray 和机器学习方面的新手。

所以我有如下 xarray 数据集:

<xarray.Dataset>
Dimensions:    (latitude: 721, longitude: 1440, time: 72)
Coordinates:
  * time       (time) datetime64[ns] 1950-01-01 1951-01-01 ... 2021-01-01
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
Data variables:
    z          (time, latitude, longitude) float32 49654.793 49654.793 ... 49654.793

现在我想跨时间维度对每个网格应用线性回归，然后我想从原始值中删除回归值以删除趋势。下面是一个示例网格的示例。

y = np.array(jan.z[:, 700, 700]) #single grid with all time
x = (np.arange(1950, len(y)+1949)).reshape(-1, 1) #72 time for x axis which will remain same for all grid

reg = LinearRegression().fit(x, y)

pred = reg.predict(x)
y = (y - (pred - y))

现在这只是一个网格，我有这样的 721*14000 网格，所以 for 循环不是最优化的方法，是否有更优化的方法或一些直接的函数所以在 xarray 中？我尝试寻找类似的东西，但找不到可以解决我的问题的东西。

Answer 1

xarray 现在有一个 polyfit 方法(参见 xr.DataArray.polyfit 和 xr.Dataset.polyfit )。它与 dask 集成良好，即使对于大型阵列也能表现得相当好。

例子

In [2]: # generate a large random dask array
   ...: arr = dda.random.random((100, 1000, 1000), chunks=(-1, 100, 100)).astype('float32')

In [3]: da = xr.DataArray(arr, dims=['time', 'y', 'x'], coords=[pd.date_range('1950-01-01', freq='Y', periods=100), np.a
  ...: range(1000), np.arange(1000)])

In [4]: da
Out[4]:
<xarray.DataArray 'astype-5e6d16cb0909caa0098a5c92f0770801' (time: 100,
                                                             y: 1000, x: 1000)>
dask.array<astype, shape=(100, 1000, 1000), dtype=float32, chunksize=(100, 100, 100), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 1950-12-31 1951-12-31 ... 2049-12-31
  * y        (y) int64 0 1 2 3 4 5 6 7 8 ... 991 992 993 994 995 996 997 998 999
  * x        (x) int64 0 1 2 3 4 5 6 7 8 ... 991 992 993 994 995 996 997 998 999

使用 deg=1 沿维度 time 调用 polyfit 进行线性最小二乘回归:

In [5]: da.polyfit('time', deg=1)
Out[5]:
<xarray.Dataset>
Dimensions:               (degree: 2, y: 1000, x: 1000)
Coordinates:
  * degree                (degree) int64 1 0
  * y                     (y) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * x                     (x) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
Data variables:
    polyfit_coefficients  (degree, y, x) float64 dask.array<chunksize=(2, 10, 1000), meta=np.ndarray>

完整的工作流程

首先，将您的年份值指定为数组的坐标，并使用 swap_dims将 'time' 替换为年度值(用于适当缩放系数)。

In [13]: year = np.arange(1950, 2050)

In [14]: da.coords['year'] = ('time', ), year

In [15]: da = da.swap_dims({'time': 'year'})

现在您可以调用 polyfit 来获取 [units/year] 的趋势:

In [16]: coeffs = da.polyfit('year', deg=1)

xr.polyval可用于使用结果系数预测值:

In [17]: pred = xr.polyval(da.year, coeffs.chunk({'x': 100, 'y': 100}))

In [18]: pred
Out[18]:
<xarray.Dataset>
Dimensions:               (y: 1000, x: 1000, year: 100)
Coordinates:
  * y                     (y) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * x                     (x) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * year                  (year) int64 1950 1951 1952 1953 ... 2047 2048 2049
Data variables:
    polyfit_coefficients  (year, y, x) float64 dask.array<chunksize=(100, 100, 100), meta=np.ndarray>

结果可用于计算无法解释的变化:

In [19]: unexplained = da - pred

In [20]: unexplained
Out[20]:
<xarray.Dataset>
Dimensions:               (y: 1000, x: 1000, year: 100)
Coordinates:
  * y                     (y) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * x                     (x) int64 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * year                  (year) int64 1950 1951 1952 1953 ... 2047 2048 2049
    time                  (year) datetime64[ns] 1950-12-31 ... 2049-12-31
Data variables:
    polyfit_coefficients  (year, y, x) float64 dask.array<chunksize=(100, 100, 100), meta=np.ndarray>