When I run the following code inside a jupyter notebook, the values on the x axis (time) in the generated plot, inserted after the code, appear to run from 0 to ~4. I expect them to run from 1 to 4, like the time values do. I can't tell if this is a problem with what xarray passes to nc_time_axis, or if it's a problem with nc_time_axis itself. Could this be looked into please?

``` import cftime import xarray as xr

time_vals = [cftime.DatetimeNoLeap(1+year, 1+month, 15) for year in range(3) for month in range(12)]

x_vals = [time_val.year + time_val.dayofyr / 365.0 for time_val in time_vals]

x_da = xr.DataArray(x_vals, coords=[time_vals], dims=["time"])

x_da.plot.line("-o"); ```

Environment:

The documentation for xarray.DataArray.copy

If deep=True, a deep copy is made of the data array. Otherwise, a shallow copy is made, so each variable in the new array’s dataset is also a variable in this array’s dataset.

I do not understand what dataset is being referred to here. In particular, there are no xarray datasets in the examples provided in this documentation. Could someone provide clarification?

Code Sample, a copy-pastable example if possible

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

create time and time_bounds DataArrays for Jan-1850 and Feb-1850

time_bounds_vals = np.array([[0.0, 31.0], [31.0, 59.0]]) time_vals = time_bounds_vals.mean(axis=1)

time_var = xr.DataArray(time_vals, dims='time', coords={'time':time_vals}) time_bounds_var = xr.DataArray(time_bounds_vals, dims=('time', 'd2'), coords={'time':time_vals})

create Dataset of time and time_bounds

ds = xr.Dataset(coords={'time':time_var}, data_vars={'time_bounds':time_bounds_var}) ds.time.attrs = {'bounds':'time_bounds', 'calendar':'noleap', 'units':'days since 1850-01-01'}

write Jan-1850 values to file

ds.isel(time=slice(0,1)).to_netcdf('Jan-1850.nc', unlimited_dims='time')

write Feb-1850 values to file

ds.isel(time=slice(1,2)).to_netcdf('Feb-1850.nc', unlimited_dims='time')

use open_mfdataset to read in files, combining into 1 Dataset

decode_times = True decode_cf = True ds = xr.open_mfdataset(['Jan-1850.nc', 'Feb-1850.nc'], decode_cf=decode_cf, decode_times=decode_times)

write combined Dataset out

ds.to_netcdf('JanFeb-1850.nc', unlimited_dims='time') ```

Problem description

The above code initially creates 2 netCDF files, for Jan-1850 and Feb-1850, that have the variables time and time_bounds, and time:bounds='time_bounds'. It then reads the 2 files back in as a single Dataset, using open_mfdataset, and this Dataset is written back out to a netCDF file. ncdump of this final file is ``` netcdf JanFeb-1850 { dimensions: time = UNLIMITED ; // (2 currently) d2 = 2 ; variables: int64 time(time) ; time:bounds = "time_bounds" ; time:units = "hours since 1850-01-16 12:00:00.000000" ; time:calendar = "noleap" ; double time_bounds(time, d2) ; time_bounds:_FillValue = NaN ; time_bounds:units = "days since 1850-01-01" ; time_bounds:calendar = "noleap" ; data:

time = 0, 708 ;

time_bounds = 0, 31, 31, 59 ; } `` The problem is that the units attribute fortimeandtime_bounds` are different in this file, contrary to what CF conventions requires.

The final call to to_netcdf is creating a file where time's units (and type) differ from what they are in the intermediate files. These transformations are not being applied to time_bounds.

While the change to time's type is not necessarily an issue, I do find it surprising.

This inconsistency goes away if either of decode_times or decode_cf is set to False in the python code above. In particular, the transformations to time's units and type do not happen.

The inconsistency also goes away if open_mfdataset opens a single file. In this scenario also, the transformations to time's units and type do not happen.

I think that the desired behavior is to either not apply the units and type transformations to time, or to also apply them to time_bounds. The first option would be consistent with the current single-file behavior.

Hi,

The code below creates a Dataset with an NxNxN DataArray that is equal to a constant val. For various re-chunked copies of the Dataset, the code computes the sum of the array, and compares it to the exact value N*N*N*val. I find that the printed values are different, at the round-off level, for different chunk sizes.

While I'm not surprised at these round-off differences, I could not find mention of such behavior in the xarray documentation.

Is this feature known to xarray developers? Do xarray developers consider it a feature or a bug?

Either way, I think it would be useful if the xarray documentation would mention that the results of some operations depends on chunk size.

code: ```import numpy as np import xarray as xr

N = 128

val = 1.9 val_array = np.full((N, N, N), val) exact_sum = N * N * N * val

ds = xr.DataArray(val_array, name='val_array', dims=['x', 'y', 'z']).to_dataset()

rel_diff = (ds['val_array'].sum().values - exact_sum) / exact_sum print('no chunking, rel_diff = %e' % rel_diff)

for chunk_x in [N//16, N//4, N]: for chunk_y in [N//16, N//4, N]: for chunk_z in [N//16, N//4, N]: ds2 = ds.chunk({'x':chunk_x, 'y':chunk_y, 'z':chunk_z}) rel_diff = (ds2['val_array'].sum().values - exact_sum) / exact_sum print('chunk_x = %3d, chunk_y = %3d, chunk_z = %3d, rel_diff = %e' \ % (chunk_x, chunk_y, chunk_z, rel_diff)) ```

results: no chunking, rel_diff = -4.557758e-15 chunk_x = 8, chunk_y = 8, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 8, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 8, chunk_z = 128, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 32, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 32, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 32, chunk_z = 128, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 128, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 128, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 8, chunk_y = 128, chunk_z = 128, rel_diff = -5.843279e-16 chunk_x = 32, chunk_y = 8, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 8, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 8, chunk_z = 128, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 32, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 32, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 32, chunk_z = 128, rel_diff = -5.843279e-16 chunk_x = 32, chunk_y = 128, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 32, chunk_y = 128, chunk_z = 32, rel_diff = -5.843279e-16 chunk_x = 32, chunk_y = 128, chunk_z = 128, rel_diff = 1.168656e-15 chunk_x = 128, chunk_y = 8, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 128, chunk_y = 8, chunk_z = 32, rel_diff = -2.337312e-16 chunk_x = 128, chunk_y = 8, chunk_z = 128, rel_diff = -5.843279e-16 chunk_x = 128, chunk_y = 32, chunk_z = 8, rel_diff = -2.337312e-16 chunk_x = 128, chunk_y = 32, chunk_z = 32, rel_diff = -5.843279e-16 chunk_x = 128, chunk_y = 32, chunk_z = 128, rel_diff = 1.168656e-15 chunk_x = 128, chunk_y = 128, chunk_z = 8, rel_diff = -5.843279e-16 chunk_x = 128, chunk_y = 128, chunk_z = 32, rel_diff = 1.168656e-15 chunk_x = 128, chunk_y = 128, chunk_z = 128, rel_diff = -4.557758e-15

Output of xr.show_versions()

It would be useful if the docs for open_dataset and open_mfdataset listed the default values of optional arguments (where there is a default).

For example, the docs for open_dataset do not list the defaults for decode_times and decode_coords, and the docs for open_mfdataset do not list the defaults for data_vars and coords.