to_zarr(append_dim="time") appends incorrect datetimes

[JackKelly]

Description

If you create a Zarr with a single timestep and then append to the time dimension of that Zarr in subsequent writes then the appended timestamps are likely to be wrong. This only seems to happen if the time dimension is datetime64.

Minimal Complete Verifiable Example

Create a really simple Dataset:

times = pd.date_range("2000-01-01 00:35", periods=8, freq="6H")
da = xr.DataArray(coords=[times], dims=["time"])
ds = da.to_dataset(name="foo")

Write just the first timestep to a new Zarr store:

ZARR_PATH = "test.zarr"
ds.isel(time=[0]).to_zarr(ZARR_PATH, mode="w")

So far, so good!

Now things get weird... let's append the remainder of ds to the Zarr store:

ds.isel(time=slice(1, None)).to_zarr(ZARR_PATH, append_dim="time")

This throws a warning, which is probably relevant:

/home/jack/miniconda3/envs/nowcasting_dataset/lib/python3.9/site-packages/xarray/core/dataset.py:2037: 
SerializationWarning: saving variable None with floating point data as an integer dtype
 without any _FillValue to use for NaNs
  return to_zarr(

What happened

Let's load the Zarr and print the contents on the time coord:

ds_loaded = xr.open_dataset(ZARR_PATH, engine="zarr")
print(ds_loaded.time)

<xarray.DataArray 'time' (time: 8)>
array(['2000-01-01T00:35', '2000-01-01T00:35',
       '2000-01-01T00:35', '2000-01-02T00:35',
       '2000-01-02T00:35', '2000-01-02T00:35',
       '2000-01-03T00:35', '2000-01-03T00:35'],
      dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2000-01-01T00:35:00 ... 2000-01-03T00:35:00

(I've removed the seconds and milliseconds to make it a bit easier to read)

The first and fifth time coords (2000-01-01T00:35 and 2000-01-02T00:35) are correct. None of the others are correct!

The encoding is not appropriate (see #3942)... notice that the units is days since..., which clearly can't represent sub-day resolution:

print(ds_loaded.time.encoding)

{'chunks': (1,),
 'preferred_chunks': {'time': 1},
 'compressor': Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0),
 'filters': None,
 'units': 'days since 2000-01-01 00:35:00',
 'calendar': 'proleptic_gregorian',
 'dtype': dtype('int64')}

What you expected to happen

The correct time coords are:

print(ds.time)

<xarray.DataArray 'time' (time: 8)>
array(['2000-01-01T00:35', '2000-01-01T06:35',
       '2000-01-01T12:35', '2000-01-01T18:35',
       '2000-01-02T00:35', '2000-01-02T06:35',
       '2000-01-02T12:35', '2000-01-02T18:35'],
      dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2000-01-01T00:35:00 ... 2000-01-02T18:35:00

Anything else we need to know?

There are three workarounds that I'm aware of:

1. When first creating the Zarr, write two or more timesteps into the Zarr. Then you can append any number of timesteps to the Zarr and everything works fine.
2. Convert the time coords to Unix epoch, represented as ints.
3. Manually set the encoding before the first write (as suggested in Time dtype encoding defaulting to int64 when writing netcdf or zarr #3942 (comment)). For example:

ds.isel(time=[0]).to_zarr(
    ZARR_PATH, 
    mode="w",
    encoding={
        'time': {
            'units': 'seconds since 1970-01-01'
        }
    }
)

Related issues

It's possible that the root cause of this issue is #3942.

And I think #3379 is another symptom of this issue.

Environment

Output of xr.show_versions()

INSTALLED VERSIONS

commit: None
python: 3.9.7 | packaged by conda-forge | (default, Sep 29 2021, 19:20:46)
[GCC 9.4.0]
python-bits: 64
OS: Linux
OS-release: 5.13.0-21-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_GB.UTF-8
LOCALE: ('en_GB', 'UTF-8')
libhdf5: 1.12.1
libnetcdf: 4.8.1

xarray: 0.20.1
pandas: 1.3.4
numpy: 1.21.4
scipy: 1.7.2
netCDF4: 1.5.8
pydap: None
h5netcdf: 0.11.0
h5py: 3.4.0
Nio: None
zarr: 2.10.1
cftime: 1.5.1.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: 1.2.8
cfgrib: 0.9.9.1
iris: None
bottleneck: 1.3.2
dask: 2021.10.0
distributed: None
matplotlib: 3.4.3
cartopy: None
seaborn: None
numbagg: None
fsspec: 2021.11.0
cupy: None
pint: None
sparse: None
setuptools: 58.5.3
pip: 21.3.1
conda: None
pytest: 6.2.5
IPython: 7.29.0
sphinx: None

[max-sixty]

Sorry this didn't get traction from a couple of years ago.

The example is high-quality! But is there a single copy-pastable example we could have to assess whether it's still a problem?

[dcherian]

Oh this is quite bad:

import pandas as pd
import xarray as xr

times = pd.date_range("2000-01-01 00:35", periods=8, freq="6h")
da = xr.DataArray(coords=[times], dims=["time"])
ds = da.to_dataset(name="foo")

ZARR_PATH = "test.zarr"
ds.isel(time=[0]).to_zarr(ZARR_PATH, mode="w")
ds.isel(time=slice(1, None)).to_zarr(ZARR_PATH, append_dim="time")
ds_loaded = xr.open_dataset(ZARR_PATH, engine="zarr")
print(ds_loaded.time)

we get

['2000-01-01T00:35:00.000000000', '2000-01-07T00:35:00.000000000',
       '2000-01-13T00:35:00.000000000', '2000-01-19T00:35:00.000000000',
       '2000-01-25T00:35:00.000000000', '2000-01-31T00:35:00.000000000',
       '2000-02-06T00:35:00.000000000', '2000-02-12T00:35:00.000000000'],
      dtype='datetime64[ns]']

instead of the expected

['2000-01-01 00:35:00', '2000-01-01 06:35:00',
               '2000-01-01 12:35:00', '2000-01-01 18:35:00',
               '2000-01-02 00:35:00', '2000-01-02 06:35:00',
               '2000-01-02 12:35:00', '2000-01-02 18:35:00'],
              dtype='datetime64[ns]'

---

We do get a warning that makes it clear what is happening though

UserWarning: Times can't be serialized faithfully to int64 with requested units 'days since 2000-01-01T00:35:00'. Serializing with units 'hours since 2000-01-01T00:35:00' instead. Set encoding['dtype'] to floating point dtype to serialize with units 'days since 2000-01-01T00:35:00'. Set encoding['units'] to 'hours since 2000-01-01T00:35:00' to silence this warning .
  ds.isel(time=slice(1, None)).to_zarr(ZARR_PATH, append_dim="time")

So we are appending values that were encoded differently

[dcherian]

Closing as dupe of #3942