Can someone explain the logic behind xarray.polyfit coefficients?

Question

I am trying to fit a linear regression to climate data from a Netcdf file. The data look like the following..

print(dsloc_lvl)

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]``

I've been using Xarray library to process data, so I've use the xarray.DataArray.polyfit and xarray.DataArray.polyval. Regression line looks good when plotting results.

However, when looking into the coefficients I've noticed they are very small. I've compare coefficients with the np.polyfit approach which are consisitent with what is expected. I figure this is because for np. ppolyfit I convert dates using date2num

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

and the xarray approach converts dates differently, I believe is with:

dsloc_lvl.time.astype(float)

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

So this makes coefficients look totally different:

np approach:

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

xarray aprroach:

dsloc_lvl.polyfit('time',1)
Out: 

Dimensions:               (degree: 2)
Coordinates:
  * degree                (degree) int32 1 0
Data variables:
    polyfit_coefficients  (degree) float64 1.525e-19 -0.1314

My question is, what are the units of time de xarray approach is using, and is there a way to scale it to match de numpy approach?

Thanks.

user_3888719 · Answer

While the results of numpy's polyfit are the regression coefficients with respect to an array of x values you pass in manually, xarray's polyfit gives coefficients in units of the coordinate labels. In the case of datetime coordinates, this often means the coefficient result is in the units of the array per nanosecond.

This happens because your data has a daily frequency but the time coordinate's labels are type datetime64[ns] (ns means nanoseconds).

Convert the linear coefficient from [1/ns] to [1/day] and you get the same result!

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

Xarray does not support numpy datetime arrays with any precision other than nanosecond, so you can't get arround this by simply changing the datetime type to, say, datetime64[D]. You can convert the coefficients you find, as above, or manually convert the axis to a float or int with the units you're looking for prior to calling polyfit.

See the xarray docs on Time Series Data for more info.

Example

As an example, I'll create a sample array:

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]:

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

If we take a look at the time coordinate, everything looks as you'd expect:

In [3]: da.time
Out[3]:

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

The trouble arises because da.polyfit needs to interpret the coordinate as numerical values. If we convert da.time to a float, you can see how we run into trouble. These values represent nanoseconds since Jan 1, 1970 0:00:00:

In [4]: da.time.astype(float)
Out[4]:

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

To get the same behavior as numpy, we could add an ordinal_day coordinate. Note here that I subtract off the start date (resulting in timedelta64[ns] data), then drop the coordinate into numpy using .values before changing precisions to timedelta64[D] (if you do this in xarray the precision change will be ignored):


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]:

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

Now we can run polyfit using ordinal_day as the coordinate (after swapping the dimensions of the array from time to ordinal_day using da.swap_dims):

In [10]: da.swap_dims({'time': 'ordinal_day'}).polyfit('ordinal_day', deg=1)
Out[10]:

Dimensions:               (degree: 2)
Coordinates:
  * degree                (degree) int64 1 0
Data variables:
    polyfit_coefficients  (degree) float64 -0.1 0.4966

This gives us the results we'd expect - I constructed the data with uniform random values in [0, 1] (so, mean 0.5 at the intercept) plus a linear trend with slope -0.1.

Can someone explain the logic behind xarray.polyfit coefficients?

1 Answers

Example