有人可以解释 xarray.polyfit 系数背后的逻辑吗？
Can someone explain the logic behind xarray.polyfit coefficients?
我正在尝试对 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 方法？
谢谢。
而 numpy's polyfit are the regression coefficients with respect to an array of x values you pass in manually, 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 之前使用您要查找的单位手动将轴转换为浮点数或整数。
有关详细信息，请参阅 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 转换为浮点数，您可以看到我们 运行 是如何陷入困境的。这些值表示自 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[D] 之前使用 .values 将坐标放入 numpy（如果你在 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
现在我们可以 运行 polyfit 使用 ordinal_day 作为坐标（在使用 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 的线性趋势构建了数据。
我正在尝试对 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 方法？
谢谢。
而 numpy's polyfit are the regression coefficients with respect to an array of x values you pass in manually, 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 之前使用您要查找的单位手动将轴转换为浮点数或整数。
有关详细信息，请参阅 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 转换为浮点数，您可以看到我们 运行 是如何陷入困境的。这些值表示自 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[D] 之前使用 .values 将坐标放入 numpy（如果你在 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
现在我们可以 运行 polyfit 使用 ordinal_day 作为坐标（在使用 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 的线性趋势构建了数据。