What should the default behaviour or polyfit be when using the time dimension

[guidocioni]

I stumbled upon this strange behaviour today and I think I've found the culprit although it still confuses me...
Suppose I have a dataset containing yearly values of days with snow on the ground considered as departure w.r.t. to the normal (1981-2010 average). I wanted to compute the linear trend in the last 30 years or so in every grid point. To start I just did some testing only using spatial averages since I noticed some weird behaviour.

If anomalies is my dataset, anomalies['snow_days'] is the variable considered with dimensions time, x, y. A quick look at the values averaged over the entire area gives me this:

anomalies.snow_days.mean(dim=['x', 'y']).plot()

Ok, makes sense, there is a clear downward trend. Let's compute that with polyfit

anomalies.snow_days.mean(dim=['x', 'y']).polyfit('time', 1, skipna=True)
> polyfit_coefficients -1.784e-17 14.42

e-17? Seems a little bit weird given the observed trend. What is going on?
Apparently the x values are converted to nanoseconds in epoch before doing the fit so the trend is actually about -0.5. I can get that if I do some transformation of the time axis

test = anomalies.assign_coords({'time': anomalies.time.dt.year - anomalies.time.dt.year[0]})
test.snow_days.mean(dim=['x', 'y']).polyfit('time', 1, skipna=True)
> polyfit_coefficients -0.563 8.232

This makes way more sense. The coefficient is now in units of days/year and not days/ns as before. However..why is the 0-order parameter so different? Looking at the plot 8 makes way more sense than 14.
Moreover, shouldn't it be stated in the documentation that, when using a time dimension for the fit, this is automatically converted to nanoseconds? I think most of the people using this function would want to compute spatial trends over a time dimension so I guess many others stumbled upon this problem.
I think the cleanest way would be to provide an additional parameter in the call to polyfit to specify the output units when using a time dimension.
What do you think? Is this expected behaviour?

[dcherian]

We could add the datetime_unit kwarg that we have for integrate and differentiate. This how it is implemented:

xarray/xarray/core/dataset.py

Lines 5945 to 5951 in 9fea799

	dim = coord_var.dims[0]
	if _contains_datetime_like_objects(coord_var):
	if coord_var.dtype.kind in "mM" and datetime_unit is None:
	datetime_unit, _ = np.datetime_data(coord_var.dtype)
	elif datetime_unit is None:
	datetime_unit = "s" # Default to seconds for cftime objects
	coord_var = coord_var._to_numeric(datetime_unit=datetime_unit)

[guidocioni]

We could add the datetime_unit kwarg that we have for integrate and differentiate. This how it is implemented:

xarray/xarray/core/dataset.py

Lines 5945 to 5951 in 9fea799

	dim = coord_var.dims[0]
	if _contains_datetime_like_objects(coord_var):
	if coord_var.dtype.kind in "mM" and datetime_unit is None:
	datetime_unit, _ = np.datetime_data(coord_var.dtype)
	elif datetime_unit is None:
	datetime_unit = "s" # Default to seconds for cftime objects
	coord_var = coord_var._to_numeric(datetime_unit=datetime_unit)

This is definitely a good idea! And the default behaviour would still be to revert to s. But you know why one of the parameter (y intercept) is completely different?