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) <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]``  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) <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  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: <xarray.Dataset> 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.

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]: <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  If we take a look at the time coordinate, everything looks as you'd expect: 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  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]: <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  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]: <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  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]: <xarray.Dataset> 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.