python - 有人可以解释 xarray.polyfit 系数背后的逻辑吗？

Question

我正在尝试对 Netcdf 文件中的气候数据进行线性回归。数据如下所示。

print(dsloc_lvl)
<xarray.DataArray 'sla' (time: 10227)>
array([0.0191, 0.0193, 0.0197, ..., 0.0936, 0.0811, 0.0695])
Coordinates:
    latitude   float32 21.62
  * time       (time) datetime64[ns] 1993-01-01 1993-01-02 ... 2020-12-31
    longitude  float32 -89.12
Attributes:
    ancillary_variables:  err_sla
    comment:              The sea level anomaly is the sea surface height abo...
    grid_mapping:         crs
    long_name:            Sea level anomaly
    standard_name:        sea_surface_height_above_sea_level
    units:                m
    _ChunkSizes:          [ 1 50 50]``

我一直在使用 Xarray 库来处理数据，所以我使用了 xarray.DataArray.polyfit 和 xarray.DataArray.polyval。绘制结果时回归线看起来不错。

但是，在查看系数时，我注意到它们非常小。我已经将系数与 np.polyfit 方法进行了比较，这与预期的一致。我认为这是因为对于 np. ppolyfit 我使用 date2num 转换日期

x1=mdates.date2num(dsloc_lvl['time'])
Out: array([ 8401.,  8402.,  8403., ..., 18625., 18626., 18627.])

并且 xarray 方法以不同方式转换日期，我相信是:

dsloc_lvl.time.astype(float)
<xarray.DataArray 'time' (time: 10227)>
array([7.2584640e+17, 7.2593280e+17, 7.2601920e+17, ..., 1.6092000e+18,
       1.6092864e+18, 1.6093728e+18])
Coordinates:
    latitude   float32 21.62
  * time       (time) datetime64[ns] 1993-01-01 1993-01-02 ... 2020-12-31
    longitude  float32 -89.12
Attributes:
    axis:                 T
    long_name:            Time
    standard_name:        time
    _ChunkSizes:          1
    _CoordinateAxisType:  Time
    valid_min:            15706.0
    valid_max:            25932.0

所以这使得系数看起来完全不同:

np 方法:

np.polyfit(x1,y1,1)
Out: array([ 1.31727420e-05, -1.31428413e-01])

xarray 方法:

dsloc_lvl.polyfit('time',1)
Out: 
<xarray.Dataset>
Dimensions:               (degree: 2)
Coordinates:
  * degree                (degree) int32 1 0
Data variables:
    polyfit_coefficients  (degree) float64 1.525e-19 -0.1314

我的问题是，de xarray 方法使用的时间单位是什么，有没有办法对其进行缩放以匹配 de numpy 方法？

谢谢。

Answer 1

虽然 numpy's polyfit 的结果是相对于您手动传入的 x 值数组的回归系数，xarray's polyfit以坐标标签为单位给出系数.在日期时间坐标的情况下，这通常意味着系数结果以数组每纳秒为单位。

发生这种情况是因为您的数据具有每日频率，但time 坐标的标签类型为datetime64[ns](ns 表示纳秒)。

将线性系数从[1/ns] 转换为[1/day]，结果相同!

1.525e-19 [units/ns] * 1e9 [ns/s] * 60 [s/m] * 60 [m/h] * 24 [h/d]
    = 1.317e-05 [units / day]

Xarray 不支持除 nanosecond 以外的任何精度的 numpy 日期时间数组，因此您无法通过简单地将日期时间类型更改为 datetime64[D]< 来解决这个问题。您可以转换找到的系数，如上所述，或者在调用 polyfit 之前使用您要查找的单位手动将轴转换为 float 或整数。

请参阅 Time Series Data 上的 xarray 文档了解更多信息。

例子

例如，我将创建一个示例数组:

In [1]: import xarray as xr, pandas as pd, numpy as np
   ...:
   ...: # create an array indexed by time, with 1096 daily observations from
   ...: # Jan 1 2020 to Dec 31, 2022. The array has noise around a linear
   ...: # trend with slope -0.1
   ...: time = pd.date_range('2020-01-01', '2022-12-31', freq='D')
   ...: Y = np.random.random(size=len(time)) + np.arange(0, (len(time) * -0.1), -0.1)
   ...: da = xr.DataArray(Y, dims=['time'], coords=[time])

In [2]: da
Out[2]:
<xarray.DataArray (time: 1096)>
array([   0.44076544,    0.66566835,    0.72999141, ..., -108.84335381,
       -109.38686183, -109.49807849])
Coordinates:
  * time     (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2022-12-31

如果我们看一下 time 坐标，一切看起来都如您所料:

In [3]: da.time
Out[3]:
<xarray.DataArray 'time' (time: 1096)>
array(['2020-01-01T00:00:00.000000000', '2020-01-02T00:00:00.000000000',
       '2020-01-03T00:00:00.000000000', ..., '2022-12-29T00:00:00.000000000',
       '2022-12-30T00:00:00.000000000', '2022-12-31T00:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2022-12-31

问题的出现是因为 da.polyfit 需要将坐标解释为数值。如果我们将 da.time 转换为 float ，您会看到我们遇到了麻烦。这些值表示自 1970 年 1 月 1 日 0:00:00 以来的纳秒数:

In [4]: da.time.astype(float)
Out[4]:
<xarray.DataArray 'time' (time: 1096)>
array([1.5778368e+18, 1.5779232e+18, 1.5780096e+18, ..., 1.6722720e+18,
       1.6723584e+18, 1.6724448e+18])
Coordinates:
  * time     (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2022-12-31

要获得与 numpy 相同的行为，我们可以添加一个 ordinal_day 坐标。请注意，我减去开始日期(产生 timedelta64[ns] 数据)，然后在将精度更改为 timedelta64 之前使用 .values 将坐标放入 numpy [D](如果您在 xarray 中执行此操作，精度更改将被忽略):

In [7]: da.coords['ordinal_day'] = (
   ...:     ('time', ),
   ...:     (da.time - da.time.min()).values.astype('timedelta64[D]').astype(int)
   ...: )

In [8]: da.ordinal_day
Out[8]:
<xarray.DataArray 'ordinal_day' (time: 1096)>
array([   0,    1,    2, ..., 1093, 1094, 1095])
Coordinates:
  * time         (time) datetime64[ns] 2020-01-01 2020-01-02 ... 2022-12-31
    ordinal_day  (time) int64 0 1 2 3 4 5 6 ... 1090 1091 1092 1093 1094 1095

现在我们可以使用 ordinal_day 作为坐标运行 polyfit(在使用 da.swap_dims 将数组的维度从 time 交换为 ordinal_day 之后) ):

In [10]: da.swap_dims({'time': 'ordinal_day'}).polyfit('ordinal_day', deg=1)
Out[10]:
<xarray.Dataset>
Dimensions:               (degree: 2)
Coordinates:
  * degree                (degree) int64 1 0
Data variables:
    polyfit_coefficients  (degree) float64 -0.1 0.4966

这给了我们预期的结果 - 我构建的数据具有 [0, 1] 中的均匀随机值(因此，截距处的平均值为 0.5)加上斜率为 -0.1 的线性趋势。