What is your issue?

When opening the xarray documentation website, it defaults to dark mode as my system is using dark mode. However, much of the text is not readable as it remains black on the dark background, see screenshots for comparison of the homepage in dark vs light mode:

Is your feature request related to a problem? Please describe. I think (please correct me if I'm wrong) when using polyfit with dim='time' the units of the output slope are [data_units]/ns, but this not explained in the docstring or in the documentation on the webpage. I figured this out eventually but it could be confusing for new users.

Describe the solution you'd like Mention in the documentation and/or docstring that the units for time will be ns

Use _possibly_convert_obhects to raise an error when converting datetime64 or timedelta64 objects that are in some units other than ns to ns as per pandas requirements, resulting in dates outside the valid range for ns precision being silently changed to incorrect values.

• [x] Closes #4427
• [x] Passes isort . && black . && mypy . && flake8
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

What happened: When using xr.DataArray.assign_coords() to assign a new coordinate to the time dimension that is an array with dtype datetime64[us], after assigning, the dtype is datetime64[ns], resulting in the wrong dates, since the dates I am using are outside the valid range for the [ns] units.

What you expected to happen: Preserve dtype of array when assigning as a coordinate.

Minimal Complete Verifiable Example:

```python import numpy as np import xarray as xr import cftime

tmp = np.random.random(12)

da = xr.DataArray(tmp, dims='time')

times=list()

for mth in np.arange(1, 13): times.append(cftime.DatetimeNoLeap(1250, mth, 1))

times64 = np.array([np.datetime64(t, 'us') for t in times])

da = da.assign_coords({'time': times64}) which gives for the original array:python In [49]: times64 Out[49]: array(['1250-01-01T00:00:00.000000', '1250-02-01T00:00:00.000000', '1250-03-01T00:00:00.000000', '1250-04-01T00:00:00.000000', '1250-05-01T00:00:00.000000', '1250-06-01T00:00:00.000000', '1250-07-01T00:00:00.000000', '1250-08-01T00:00:00.000000', '1250-09-01T00:00:00.000000', '1250-10-01T00:00:00.000000', '1250-11-01T00:00:00.000000', '1250-12-01T00:00:00.000000'], dtype='datetime64[us]') ```

and for the array after assigning: python In [51]: da.time Out[51]: <xarray.DataArray 'time' (time: 12)> array(['1834-07-22T23:34:33.709551616', '1834-08-22T23:34:33.709551616', '1834-09-19T23:34:33.709551616', '1834-10-20T23:34:33.709551616', '1834-11-19T23:34:33.709551616', '1834-12-20T23:34:33.709551616', '1835-01-19T23:34:33.709551616', '1835-02-19T23:34:33.709551616', '1835-03-22T23:34:33.709551616', '1835-04-21T23:34:33.709551616', '1835-05-22T23:34:33.709551616', '1835-06-21T23:34:33.709551616'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 1834-07-22T23:34:33.709551616 ... 1835-06-...

Anything else we need to know?:

Environment:

MCVE Code Sample

```python import numpy as np import xarray as xr import cftime

create time coordinate

tdays = np.arange(0, 730) time = cftime.num2date(tdays, 'days since 0001-01-01 00:00:00', calendar='noleap')

create spatial coordinate

lev = np.arange(100)

Create dummy data

x = np.random.rand(time.size, lev.size) y = np.random.rand(lev.size)

Create sample Dataset

ds = xr.Dataset({'sample_data': (['time', 'lev'], x), 'independent_data': (['lev'], y)}, coords={'time': (['time'], time), 'lev': (['lev'], lev)})

Perform groupby and mean

ds2 = ds.groupby('time.month').mean(dim='time') ```

Actual Output

python <xarray.Dataset> Dimensions: (lev: 100, month: 12) Coordinates: * lev (lev) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99 * month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12 Data variables: sample_data (month, lev) float64 0.5143 0.554 0.5027 ... 0.5246 0.5435 independent_data (month, lev) float64 0.01667 0.4687 ... 0.1015 0.7459

Expected Output

python <xarray.Dataset> Dimensions: (lev: 100, month: 12) Coordinates: * lev (lev) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99 * month (month) int64 1 2 3 4 5 6 7 8 9 10 11 12 Data variables: sample_data (month, lev) float64 0.5143 0.554 0.5027 ... 0.5246 0.5435 independent_data (lev) float64 0.01667 0.4687 ... 0.1015 0.7459

Problem Description

The variable independent_data above initially has no time dimension but, after performing groupby('time.month').mean(dim='time') on the Dataset, it now has a month dimension that is meaningless. Preferably, it should leave the independent_data variable untouched.

Output of xr.show_versions()