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?

MCVE Code Sample

``` import numpy as np import xarray as xr

xlen = 4 x = xr.DataArray(np.linspace(0.0, 1.0, xlen), dims=('x'))

varname = 'foo' xr.Dataset({varname: xr.DataArray(np.arange(xlen, dtype='float64'), dims=('x'), coords={'x': x})}).to_netcdf('ds.nc')

with xr.open_dataset('ds.nc') as ds: ds2 = ds.copy() ds2[varname][0] = 11.0 print(f'ds.equals = {ds.equals(ds2)}')

with xr.open_dataset('ds.nc') as ds: ds2 = ds.copy() print(f'ds.equals = {ds.equals(ds2)}') ds2[varname][0] = 11.0 print(f'ds.equals = {ds.equals(ds2)}') ```

Expected Output

I expect the code to write out ds.equals = True ds.equals = True ds.equals = True

However, when I run it, the last line is ds.equals = False

Problem Description

The code above writes a small xr.Dataset to a netCDF file. There are 2 context managers opening the netCDF file as ds. Both context manager blocks start by setting ds2 to a shallow copy of ds.

In the first context manager block, a value in ds2 is modified, and ds2 is compared to ds. The Datasets are still equal, confirming that the copy is shallow.

The second context manager block is the same as the first, except that ds2 is compared to ds prior changing the value the value in ds2. When this is done, the Datasets are no longer equal, indicating that ds2 is no longer a shallow copy of ds.

I don't understand how evaluating ds.equals(ds2), prior to changing a value in ds2, could decouple ds2 from ds.

I only observe this behavior when ds is set via xr.open_dataset. I don't see it when I create ds directly using xr.Dataset.

I'm rather perplexed by this.

Output of xr.show_versions()

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.