Repository: pydata/xarray

Issue:
Differences on datetime values appears after writing reindexed variable on netCDF file
In my Dataset i've got a time serie coordinate who begins like this 

```
<xarray.DataArray 'time' (time: 10)>
array(['2014-02-15T00:00:00.000000000+0100',
       '2014-02-15T18:10:00.000000000+0100',
       '2014-02-16T18:10:00.000000000+0100',
       '2014-02-17T18:10:00.000000000+0100',
       '2014-02-18T18:10:00.000000000+0100',
       '2014-02-19T18:10:00.000000000+0100',
       '2014-02-20T18:10:00.000000000+0100',
       '2014-02-21T18:10:00.000000000+0100',
       '2014-02-22T00:00:00.000000000+0100',
       '2014-02-23T00:00:00.000000000+0100'], dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2014-02-14T23:00:00 2014-02-15T17:10:00 ...
```

And all is ok when I write and re-open the netdcdf file

Then i try to add to this dataset a reindexed variable like this

`da["MeanRainfallHeigh"] = rain.reindex(time =da.time).fillna(0)`

Everything is still good for the writing, but when I reopen the netcdf file, the time values are modified for the minutes part.

```
<xarray.DataArray 'time' (time: 10)>
array(['2014-02-15T00:00:00.000000000+0100',
       '2014-02-15T18:00:00.000000000+0100',
       '2014-02-16T18:00:00.000000000+0100',
       '2014-02-17T18:00:00.000000000+0100',
       '2014-02-18T18:00:00.000000000+0100',
       '2014-02-19T18:00:00.000000000+0100',
       '2014-02-20T18:00:00.000000000+0100',
       '2014-02-21T18:00:00.000000000+0100',
       '2014-02-22T00:00:00.000000000+0100',
       '2014-02-23T00:00:00.000000000+0100'], dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2014-02-14T23:00:00 2014-02-15T17:00:00 ...
```

Thanks!


Discussion:
kmuehlbauer: Unfortunately this is more involved. Coming back later.
headtr1ck: > Unfortunately this is more involved. Coming back later.

Ah sry, didn't see this. Good luck
kmuehlbauer: > > Unfortunately this is more involved. Coming back later.
> 
> Ah sry, didn't see this. Good luck

Thanks, finally, that wasn't that complicated :-)
dcherian: Should we override here or warn the user instead?
kmuehlbauer: > Should we override here or warn the user instead?

Quoting from @NotSqrt's comment https://github.com/pydata/xarray/issues/1064#issuecomment-1723248713:

_If the resolution loss can't be fixed automatically, what would be nice in the warning is a link or a summary of what the user has to do to solve the resolution loss!_

I'm totally open here. The question is, what takes precedence `dtype` or `units`? It doesn't make much sense to give `units` which contradict `dtype` (at least for faithful write). So these kind of issues mostly (that's debatable) are not because the user giving problematic encoding parameters, but due to inconsistencies while processing. 
  
kmuehlbauer: To add to my above comment, with #7827 we added a warning when times/timedeltas had to be encoded to float64 because of units not fitting.

In this PR for those cases where times/timedeltas had to be encoded to fit `units` and `encoding["dtype"] != np.float64` and of integer dtype `encoding["dtype"] will be dropped to make further cf encoding work. 

If needed we could add information to the above mentioned warning that `encoding["dtype"]` will be overridden or that dtype/units combination does not fit. For that case it might be necessary to move that warning out of `encode_cf_datetime`.

Which approach should we take?



spencerkclark: Yeah, it seems like the two options are:
- Warn and override (to @dcherian's point, it probably would be useful to include information that the encoding dtype is being overridden in the warning).
- Raise an informative error if the specified encoding units and dtype are incompatible with each other.

> The question is, what takes precedence dtype or units? It doesn't make much sense to give units which contradict dtype (at least for faithful write).

Exactly, I think this is the right way to frame it.  I guess the answer is really up to the user.  It's hard to know ahead of time what they care about more, so maybe the most conservative approach is to raise?  Or do we feel having a default solution (with information on how to address it in other ways) is more helpful?

NotSqrt: To give a little context, I initially encountered those problems in this setup:

- I receive some data (similar to MQTT), and aggregate it in an xarray + serialized as netcdf.
- the time between each data can be quite random, but it may happen that the first points of the file are separated by a round number of minutes. Which means that the serialization will decide to use 'minutes since REF_TIME'
- when the subsequent data point is no longer separated by a round number of minutes, but only some seconds, bam, bug, the seconds were lost when aggregated into an existing file.

I don't care about the exact unit, I just don't want to randomly lose some precision on my data ..
kmuehlbauer: > It's hard to know ahead of time what they care about more, so maybe the most conservative approach is to raise? Or do we feel having a default solution (with information on how to address it in other ways) is more helpful?

@spencerkclark @dcherian Yes, it's probably best to raise with a meaningful error message. 



NotSqrt: > > It's hard to know ahead of time what they care about more, so maybe the most conservative approach is to raise? Or do we feel having a default solution (with information on how to address it in other ways) is more helpful?
> 
> @spencerkclark @dcherian Yes, it's probably best to raise with a meaningful error message.

Having a random exception pop up in production just because the delta between points has changed would have been very surprising.
At the time, I expected the serialization to automatically handle that it had to switch from "minutes since REF_TIME" into "seconds since REF_TIME"
dcherian: There is no way to decide if `encoding` was set in a file and not modified later. 

>  I guess the answer is really up to the user. It's hard to know ahead of time what they care about more, so maybe the most conservative approach is to raise?

Yes, I agree with this framing. 

kmuehlbauer: @dcherian @spencerkclark I'm a bit biased here because my own workflows have similar issues as @NotSqrt's.

I'm now thinking that my first approach overriding `dtype` did not work out. 

If provided `units` would require `dtype` change for faithful serialization why not just fix `units`, instead of raising? We will warn the user that `units` will be changed and can keep `dtype`.

kmuehlbauer: The last iteration aligns with my case 1 from above.

This is the flow now:

1. check if times can be serialized faithfully with given `units` into int64.
yes? -> all good, serialize to int64
no? -> follow 2
2. check encoding dtype:
  a. None -> `units` take precedence, just warn, serialize to float64
  b. integer -> warn and override `units` to preserve integer `dtype` (and with that nanosecond resolution) -> serialize to int64
  c. else (floating) -> serialize to float64 without warning
  
We could raise instead of warn in 2b, but I'm not sure that will do any good for automated workflows.


spencerkclark: Thanks @kmuehlbauer.  I agree, if we are going to override by default, overriding the units is more palatable than the dtype, since it guarantees round-trip accuracy (in other words if one is staying in the xarray ecosystem, it's fairly clearly the way to go).  

While I can come up with hypothetical reasons for why one might want to preserve units encoding over the dtype, I don't think it would be a very common situation.  So I'd say I am on board with your current approach, though we should make sure that the warning describes how one can preserve the units instead of the dtype by modifying the encoding if they'd like a different behavior.
kmuehlbauer: @spencerkclark I've added the needed instructions for achieving different behaviour, plus howto silence these warning. This should be good to go then from my side.

It would be good to get this into the next version alongside the nanosecond fix. 

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 67429ed7..5f18e999 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -90,6 +90,8 @@ Bug fixes
 - ``.rolling_exp`` functions no longer mistakenly lose non-dimensioned coords
   (:issue:`6528`, :pull:`8114`)
   By `Maximilian Roos <https://github.com/max-sixty>`_.
+- In the event that user-provided datetime64/timedelta64 units and integer dtype encoding parameters conflict with each other, override the units to preserve an integer dtype for most faithful serialization to disk (:issue:`1064`, :pull:`8201`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index 79efbecf..2822f02d 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -656,8 +656,22 @@ def cast_to_int_if_safe(num) -> np.ndarray:
     return num
 
 
+def _division(deltas, delta, floor):
+    if floor:
+        # calculate int64 floor division
+        # to preserve integer dtype if possible (GH 4045, GH7817).
+        num = deltas // delta.astype(np.int64)
+        num = num.astype(np.int64, copy=False)
+    else:
+        num = deltas / delta
+    return num
+
+
 def encode_cf_datetime(
-    dates, units: str | None = None, calendar: str | None = None
+    dates,
+    units: str | None = None,
+    calendar: str | None = None,
+    dtype: np.dtype | None = None,
 ) -> tuple[np.ndarray, str, str]:
     """Given an array of datetime objects, returns the tuple `(num, units,
     calendar)` suitable for a CF compliant time variable.
@@ -689,6 +703,12 @@ def encode_cf_datetime(
         time_units, ref_date = _unpack_time_units_and_ref_date(units)
         time_delta = _time_units_to_timedelta64(time_units)
 
+        # Wrap the dates in a DatetimeIndex to do the subtraction to ensure
+        # an OverflowError is raised if the ref_date is too far away from
+        # dates to be encoded (GH 2272).
+        dates_as_index = pd.DatetimeIndex(dates.ravel())
+        time_deltas = dates_as_index - ref_date
+
         # retrieve needed units to faithfully encode to int64
         needed_units, data_ref_date = _unpack_time_units_and_ref_date(data_units)
         if data_units != units:
@@ -697,26 +717,32 @@ def encode_cf_datetime(
             if ref_delta > np.timedelta64(0, "ns"):
                 needed_units = _infer_time_units_from_diff(ref_delta)
 
-        # Wrap the dates in a DatetimeIndex to do the subtraction to ensure
-        # an OverflowError is raised if the ref_date is too far away from
-        # dates to be encoded (GH 2272).
-        dates_as_index = pd.DatetimeIndex(dates.ravel())
-        time_deltas = dates_as_index - ref_date
-
         # needed time delta to encode faithfully to int64
         needed_time_delta = _time_units_to_timedelta64(needed_units)
-        if time_delta <= needed_time_delta:
-            # calculate int64 floor division
-            # to preserve integer dtype if possible (GH 4045, GH7817).
-            num = time_deltas // time_delta.astype(np.int64)
-            num = num.astype(np.int64, copy=False)
-        else:
-            emit_user_level_warning(
-                f"Times can't be serialized faithfully with requested units {units!r}. "
-                f"Resolution of {needed_units!r} needed. "
-                f"Serializing timeseries to floating point."
-            )
-            num = time_deltas / time_delta
+
+        floor_division = True
+        if time_delta > needed_time_delta:
+            floor_division = False
+            if dtype is None:
+                emit_user_level_warning(
+                    f"Times can't be serialized faithfully to int64 with requested units {units!r}. "
+                    f"Resolution of {needed_units!r} needed. Serializing times to floating point instead. "
+                    f"Set encoding['dtype'] to integer dtype to serialize to int64. "
+                    f"Set encoding['dtype'] to floating point dtype to silence this warning."
+                )
+            elif np.issubdtype(dtype, np.integer):
+                new_units = f"{needed_units} since {format_timestamp(ref_date)}"
+                emit_user_level_warning(
+                    f"Times can't be serialized faithfully to int64 with requested units {units!r}. "
+                    f"Serializing with units {new_units!r} instead. "
+                    f"Set encoding['dtype'] to floating point dtype to serialize with units {units!r}. "
+                    f"Set encoding['units'] to {new_units!r} to silence this warning ."
+                )
+                units = new_units
+                time_delta = needed_time_delta
+                floor_division = True
+
+        num = _division(time_deltas, time_delta, floor_division)
         num = num.values.reshape(dates.shape)
 
     except (OutOfBoundsDatetime, OverflowError, ValueError):
@@ -728,7 +754,9 @@ def encode_cf_datetime(
     return (num, units, calendar)
 
 
-def encode_cf_timedelta(timedeltas, units: str | None = None) -> tuple[np.ndarray, str]:
+def encode_cf_timedelta(
+    timedeltas, units: str | None = None, dtype: np.dtype | None = None
+) -> tuple[np.ndarray, str]:
     data_units = infer_timedelta_units(timedeltas)
 
     if units is None:
@@ -744,18 +772,29 @@ def encode_cf_timedelta(timedeltas, units: str | None = None) -> tuple[np.ndarra
 
     # needed time delta to encode faithfully to int64
     needed_time_delta = _time_units_to_timedelta64(needed_units)
-    if time_delta <= needed_time_delta:
-        # calculate int64 floor division
-        # to preserve integer dtype if possible
-        num = time_deltas // time_delta.astype(np.int64)
-        num = num.astype(np.int64, copy=False)
-    else:
-        emit_user_level_warning(
-            f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
-            f"Resolution of {needed_units!r} needed. "
-            f"Serializing timedeltas to floating point."
-        )
-        num = time_deltas / time_delta
+
+    floor_division = True
+    if time_delta > needed_time_delta:
+        floor_division = False
+        if dtype is None:
+            emit_user_level_warning(
+                f"Timedeltas can't be serialized faithfully to int64 with requested units {units!r}. "
+                f"Resolution of {needed_units!r} needed. Serializing timeseries to floating point instead. "
+                f"Set encoding['dtype'] to integer dtype to serialize to int64. "
+                f"Set encoding['dtype'] to floating point dtype to silence this warning."
+            )
+        elif np.issubdtype(dtype, np.integer):
+            emit_user_level_warning(
+                f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
+                f"Serializing with units {needed_units!r} instead. "
+                f"Set encoding['dtype'] to floating point dtype to serialize with units {units!r}. "
+                f"Set encoding['units'] to {needed_units!r} to silence this warning ."
+            )
+            units = needed_units
+            time_delta = needed_time_delta
+            floor_division = True
+
+    num = _division(time_deltas, time_delta, floor_division)
     num = num.values.reshape(timedeltas.shape)
     return (num, units)
 
@@ -772,7 +811,8 @@ class CFDatetimeCoder(VariableCoder):
 
             units = encoding.pop("units", None)
             calendar = encoding.pop("calendar", None)
-            (data, units, calendar) = encode_cf_datetime(data, units, calendar)
+            dtype = encoding.get("dtype", None)
+            (data, units, calendar) = encode_cf_datetime(data, units, calendar, dtype)
 
             safe_setitem(attrs, "units", units, name=name)
             safe_setitem(attrs, "calendar", calendar, name=name)
@@ -807,7 +847,9 @@ class CFTimedeltaCoder(VariableCoder):
         if np.issubdtype(variable.data.dtype, np.timedelta64):
             dims, data, attrs, encoding = unpack_for_encoding(variable)
 
-            data, units = encode_cf_timedelta(data, encoding.pop("units", None))
+            data, units = encode_cf_timedelta(
+                data, encoding.pop("units", None), encoding.get("dtype", None)
+            )
             safe_setitem(attrs, "units", units, name=name)
 
             return Variable(dims, data, attrs, encoding, fastpath=True)


Test patch:
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index 079e432b..5f76a4a2 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -30,7 +30,7 @@ from xarray.coding.times import (
 from xarray.coding.variables import SerializationWarning
 from xarray.conventions import _update_bounds_attributes, cf_encoder
 from xarray.core.common import contains_cftime_datetimes
-from xarray.testing import assert_allclose, assert_equal, assert_identical
+from xarray.testing import assert_equal, assert_identical
 from xarray.tests import (
     FirstElementAccessibleArray,
     arm_xfail,
@@ -1036,7 +1036,7 @@ def test_encode_cf_datetime_defaults_to_correct_dtype(
         pytest.skip("Nanosecond frequency is not valid for cftime dates.")
     times = date_range("2000", periods=3, freq=freq)
     units = f"{encoding_units} since 2000-01-01"
-    encoded, _, _ = coding.times.encode_cf_datetime(times, units)
+    encoded, _units, _ = coding.times.encode_cf_datetime(times, units)
 
     numpy_timeunit = coding.times._netcdf_to_numpy_timeunit(encoding_units)
     encoding_units_as_timedelta = np.timedelta64(1, numpy_timeunit)
@@ -1212,6 +1212,7 @@ def test_contains_cftime_lazy() -> None:
         ("1677-09-21T00:12:43.145224193", "ns", np.int64, None, False),
         ("1677-09-21T00:12:43.145225", "us", np.int64, None, False),
         ("1970-01-01T00:00:01.000001", "us", np.int64, None, False),
+        ("1677-09-21T00:21:52.901038080", "ns", np.float32, 20.0, True),
     ],
 )
 def test_roundtrip_datetime64_nanosecond_precision(
@@ -1261,14 +1262,52 @@ def test_roundtrip_datetime64_nanosecond_precision_warning() -> None:
     ]
     units = "days since 1970-01-10T01:01:00"
     needed_units = "hours"
-    encoding = dict(_FillValue=20, units=units)
+    new_units = f"{needed_units} since 1970-01-10T01:01:00"
+
+    encoding = dict(dtype=None, _FillValue=20, units=units)
     var = Variable(["time"], times, encoding=encoding)
-    wmsg = (
-        f"Times can't be serialized faithfully with requested units {units!r}. "
-        f"Resolution of {needed_units!r} needed. "
-    )
-    with pytest.warns(UserWarning, match=wmsg):
+    with pytest.warns(UserWarning, match=f"Resolution of {needed_units!r} needed."):
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.float64
+    assert encoded_var.attrs["units"] == units
+    assert encoded_var.attrs["_FillValue"] == 20.0
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="int64", _FillValue=20, units=units)
+    var = Variable(["time"], times, encoding=encoding)
+    with pytest.warns(
+        UserWarning, match=f"Serializing with units {new_units!r} instead."
+    ):
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == new_units
+    assert encoded_var.attrs["_FillValue"] == 20
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="float64", _FillValue=20, units=units)
+    var = Variable(["time"], times, encoding=encoding)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.float64
+    assert encoded_var.attrs["units"] == units
+    assert encoded_var.attrs["_FillValue"] == 20.0
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+    encoding = dict(dtype="int64", _FillValue=20, units=new_units)
+    var = Variable(["time"], times, encoding=encoding)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
         encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == new_units
+    assert encoded_var.attrs["_FillValue"] == 20
 
     decoded_var = conventions.decode_cf_variable("foo", encoded_var)
     assert_identical(var, decoded_var)
@@ -1309,14 +1348,18 @@ def test_roundtrip_timedelta64_nanosecond_precision_warning() -> None:
     needed_units = "hours"
     wmsg = (
         f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
-        f"Resolution of {needed_units!r} needed. "
+        f"Serializing with units {needed_units!r} instead."
     )
-    encoding = dict(_FillValue=20, units=units)
+    encoding = dict(dtype=np.int64, _FillValue=20, units=units)
     var = Variable(["time"], timedelta_values, encoding=encoding)
     with pytest.warns(UserWarning, match=wmsg):
         encoded_var = conventions.encode_cf_variable(var)
+    assert encoded_var.dtype == np.int64
+    assert encoded_var.attrs["units"] == needed_units
+    assert encoded_var.attrs["_FillValue"] == 20
     decoded_var = conventions.decode_cf_variable("foo", encoded_var)
-    assert_allclose(var, decoded_var)
+    assert_identical(var, decoded_var)
+    assert decoded_var.encoding["dtype"] == np.int64
 
 
 def test_roundtrip_float_times() -> None: