Repository: pydata/xarray

Issue:
writing datetime64 in netCDF may produce badly formatted or unreadable files
### What happened?

When writing datetime64 and NaT values  to netCDF (using `h5netcdf` engine) using user defined dtype and _FillValue (to comply to some format requirements), I get different results depending on the precision of datetime64 values.

Here I'm writing datetime64 with milliseconds precision:

```
import xarray as xr
import numpy as np

# writing datetime64 with milliseconds
times = np.ma.array([
   np.datetime64('2010-01-01 12:00:00.005'), np.datetime64('NaT')]).astype('datetime64[ns]')

da = xr.DataArray(
    name='time',
    data=times,
    dims=('time',))
da.encoding['units'] = "seconds since 2010-01-01 00:00:00.000"
da.encoding['_FillValue'] = 1e20
da.encoding['dtype'] = np.float64
ds=xr.Dataset(coords={'time': da})
ds.to_netcdf('milliseconds.nc', engine='h5netcdf')

new_ds = xr.open_dataset('milliseconds.nc')
print(new_ds.time)
```
above works fine, and fill values are correctly written, as shown by the CDL output:

```
netcdf milliseconds {
dimensions:
	time = 2 ;
variables:
	double time(time) ;
		time:_FillValue = 1.e+20 ;
		string time:units = "seconds since 2010-01-01" ;
		string time:calendar = "proleptic_gregorian" ;
data:

 time = 43200.005, _ ;
}
```

However, if instead I write datetime64 rounded to seconds, the output result is wrong:

```
# writing datetime64 with seconds
times = np.ma.array([
    np.datetime64('2010-01-01 12:00:00'), np.datetime64('NaT')]).astype('datetime64[ns]')

da = xr.DataArray(
    name='time',
    data=times,
    dims=('time',))
da.encoding['units'] = "seconds since 2010-01-01 00:00:00.000"
da.encoding['_FillValue'] = 1e20
da.encoding['dtype'] = np.float64
ds=xr.Dataset(coords={'time': da})
ds.to_netcdf('seconds.nc', engine='h5netcdf')

new_ds = xr.open_dataset('seconds.nc')
```

The CDL output is as follow (note that -9.22337203685478e+18 is used instead of the required 1e20 as _FillValue):

```
netcdf seconds {
dimensions:
	time = 2 ;
variables:
	double time(time) ;
		time:_FillValue = 1.e+20 ;
		string time:units = "seconds since 2010-01-01" ;
		string time:calendar = "proleptic_gregorian" ;
data:

 time = 43200, -9.22337203685478e+18 ;
}
```

As the result, attempting to read the result file will fail:
```
Traceback (most recent call last):
  File "conversion.pyx", line 228, in pandas._libs.tslibs.conversion.cast_from_unit
OverflowError: Python int too large to convert to C long

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "timedeltas.pyx", line 377, in pandas._libs.tslibs.timedeltas._maybe_cast_from_unit
  File "conversion.pyx", line 230, in pandas._libs.tslibs.conversion.cast_from_unit
pandas._libs.tslibs.np_datetime.OutOfBoundsDatetime: cannot convert input -9.223372036854776e+18 with the unit 's'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 326, in decode_cf_datetime
    dates = _decode_datetime_with_pandas(flat_num_dates, units, calendar)
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 272, in _decode_datetime_with_pandas
    pd.to_timedelta(flat_num_dates.min(), time_units) + ref_date
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/pandas/core/tools/timedeltas.py", line 223, in to_timedelta
    return _coerce_scalar_to_timedelta_type(arg, unit=unit, errors=errors)
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/pandas/core/tools/timedeltas.py", line 233, in _coerce_scalar_to_timedelta_type
    result = Timedelta(r, unit)
  File "timedeltas.pyx", line 1896, in pandas._libs.tslibs.timedeltas.Timedelta.__new__
  File "timedeltas.pyx", line 356, in pandas._libs.tslibs.timedeltas.convert_to_timedelta64
  File "timedeltas.pyx", line 379, in pandas._libs.tslibs.timedeltas._maybe_cast_from_unit
pandas._libs.tslibs.np_datetime.OutOfBoundsTimedelta: Cannot cast -9.223372036854776e+18 from s to 'ns' without overflow.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 218, in _decode_cf_datetime_dtype
    result = decode_cf_datetime(example_value, units, calendar, use_cftime)
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 328, in decode_cf_datetime
    dates = _decode_datetime_with_cftime(
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 242, in _decode_datetime_with_cftime
    cftime.num2date(num_dates, units, calendar, only_use_cftime_datetimes=True)
  File "src/cftime/_cftime.pyx", line 627, in cftime._cftime.num2date
  File "src/cftime/_cftime.pyx", line 421, in cftime._cftime.cast_to_int
OverflowError: time values outside range of 64 bit signed integers

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/conventions.py", line 450, in decode_cf_variables
    new_vars[k] = decode_cf_variable(
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/conventions.py", line 291, in decode_cf_variable
    var = times.CFDatetimeCoder(use_cftime=use_cftime).decode(var, name=name)
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 992, in decode
    dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/coding/times.py", line 228, in _decode_cf_datetime_dtype
    raise ValueError(msg)
ValueError: unable to decode time units 'seconds since 2010-01-01' with "calendar 'proleptic_gregorian'". Try opening your dataset with decode_times=False or installing cftime if it is not installed.

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/jfpiolle/git/felyx_processor/test_time.py", line 39, in <module>
    new_ds = xr.open_dataset('seconds.nc')
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/backends/api.py", line 588, in open_dataset
    backend_ds = backend.open_dataset(
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/backends/h5netcdf_.py", line 418, in open_dataset
    ds = store_entrypoint.open_dataset(
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/backends/store.py", line 46, in open_dataset
    vars, attrs, coord_names = conventions.decode_cf_variables(
  File "/home/jfpiolle/mambaforge/envs/poetry/lib/python3.10/site-packages/xarray/conventions.py", line 461, in decode_cf_variables
    raise type(e)(f"Failed to decode variable {k!r}: {e}") from e
ValueError: Failed to decode variable 'time': unable to decode time units 'seconds since 2010-01-01' with "calendar 'proleptic_gregorian'". Try opening your dataset with decode_times=False or installing cftime if it is not installed.
```

### What did you expect to happen?

the content of the output in the second example (with values rounded to seconds) should be similar to the first one (minus the precision of the values of course) since same units, _FillValue and variable dtype are used.

I did not met this issue with earlier versions of xarray.

### Minimal Complete Verifiable Example

```Python
see in issue description.
```


### MVCE confirmation

- [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
- [ ] Complete example — the example is self-contained, including all data and the text of any traceback.
- [ ] Verifiable example — the example copy & pastes into an IPython prompt or [Binder notebook](https://mybinder.org/v2/gh/pydata/xarray/main?urlpath=lab/tree/doc/examples/blank_template.ipynb), returning the result.
- [ ] New issue — a search of GitHub Issues suggests this is not a duplicate.
- [ ] Recent environment — the issue occurs with the latest version of xarray and its dependencies.

### Relevant log output

_No response_

### Anything else we need to know?

_No response_

### Environment

Using:
```
python                    3.10.14         hd12c33a_0_cpython    conda-forge
xarray                    2024.7.0                 pypi_0    pypi
numpy                     2.1.1                    pypi_0    pypi
h5netcdf                  1.3.0                    pypi_0    pypi
h5py                      3.11.0                   pypi_0    pypi
```

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
python-bits: 64
OS: Linux
OS-release: 5.4.0-150-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: fr_FR.UTF-8
LOCALE: ('fr_FR', 'UTF-8')
libhdf5: 1.14.2
libnetcdf: None

xarray: 2024.7.0
pandas: 2.2.2
numpy: 2.1.1
scipy: 1.14.1
netCDF4: None
pydap: None
h5netcdf: 1.3.0
h5py: 3.11.0
zarr: None
cftime: 1.6.4
nc_time_axis: None
iris: None
bottleneck: None
dask: 2024.8.2
distributed: None
matplotlib: None
cartopy: None
seaborn: None
numbagg: None
fsspec: 2024.9.0
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 74.1.2
pip: 24.2
conda: None
pytest: 6.2.5
mypy: None
IPython: None
sphinx: None

</details>


Discussion:
kmuehlbauer: Awesome, really just right around the corner. Thanks so much @spencerkclark!

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index e6caac78..264c07f5 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -30,7 +30,8 @@ New Features
   `Matt Savoie <https://github.com/flamingbear>`_,
   `Stephan Hoyer <https://github.com/shoyer>`_ and
   `Tom Nicholas <https://github.com/TomNicholas>`_.
-
+- Added zarr backends for :py:func:`open_groups` (:issue:`9430`, :pull:`9469`).
+  By `Eni Awowale <https://github.com/eni-awowale>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
@@ -46,7 +47,10 @@ Bug fixes
 - Make illegal path-like variable names when constructing a DataTree from a Dataset
   (:issue:`9339`, :pull:`9378`)
   By `Etienne Schalk <https://github.com/etienneschalk>`_.
-
+- Fix bug when encoding times with missing values as floats in the case when
+  the non-missing times could in theory be encoded with integers
+  (:issue:`9488`, :pull:`9497`). By `Spencer Clark
+  <https://github.com/spencerkclark>`_.
 
 
 Documentation
@@ -182,8 +186,8 @@ New Features
   to return an object without ``attrs``. A ``deep`` parameter controls whether
   variables' ``attrs`` are also dropped.
   By `Maximilian Roos <https://github.com/max-sixty>`_. (:pull:`8288`)
+- Added :py:func:`open_groups` for h5netcdf and netCDF4 backends (:issue:`9137`, :pull:`9243`).
   By `Eni Awowale <https://github.com/eni-awowale>`_.
-- Add `open_groups` method for unaligned datasets (:issue:`9137`, :pull:`9243`)
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index cfdecd28..5655bd20 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -771,7 +771,7 @@ def _eagerly_encode_cf_datetime(
         # needed time delta to encode faithfully to int64
         needed_time_delta = _time_units_to_timedelta64(needed_units)
 
-        floor_division = True
+        floor_division = np.issubdtype(dtype, np.integer) or dtype is None
         if time_delta > needed_time_delta:
             floor_division = False
             if dtype is None:
@@ -892,7 +892,7 @@ def _eagerly_encode_cf_timedelta(
     # needed time delta to encode faithfully to int64
     needed_time_delta = _time_units_to_timedelta64(needed_units)
 
-    floor_division = True
+    floor_division = np.issubdtype(dtype, np.integer) or dtype is None
     if time_delta > needed_time_delta:
         floor_division = False
         if dtype is None:


Test patch:
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index 5879e6be..66caf25c 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -1383,16 +1383,38 @@ def test_roundtrip_timedelta64_nanosecond_precision_warning() -> None:
     assert decoded_var.encoding["dtype"] == np.int64
 
 
-def test_roundtrip_float_times() -> None:
-    # Regression test for GitHub issue #8271
-    fill_value = 20.0
-    times = [
-        np.datetime64("1970-01-01 00:00:00", "ns"),
-        np.datetime64("1970-01-01 06:00:00", "ns"),
-        np.datetime64("NaT", "ns"),
-    ]
+_TEST_ROUNDTRIP_FLOAT_TIMES_TESTS = {
+    "GH-8271": (
+        20.0,
+        np.array(
+            ["1970-01-01 00:00:00", "1970-01-01 06:00:00", "NaT"],
+            dtype="datetime64[ns]",
+        ),
+        "days since 1960-01-01",
+        np.array([3653, 3653.25, 20.0]),
+    ),
+    "GH-9488-datetime64[ns]": (
+        1.0e20,
+        np.array(["2010-01-01 12:00:00", "NaT"], dtype="datetime64[ns]"),
+        "seconds since 2010-01-01",
+        np.array([43200, 1.0e20]),
+    ),
+    "GH-9488-timedelta64[ns]": (
+        1.0e20,
+        np.array([1_000_000_000, "NaT"], dtype="timedelta64[ns]"),
+        "seconds",
+        np.array([1.0, 1.0e20]),
+    ),
+}
+
 
-    units = "days since 1960-01-01"
+@pytest.mark.parametrize(
+    ("fill_value", "times", "units", "encoded_values"),
+    _TEST_ROUNDTRIP_FLOAT_TIMES_TESTS.values(),
+    ids=_TEST_ROUNDTRIP_FLOAT_TIMES_TESTS.keys(),
+)
+def test_roundtrip_float_times(fill_value, times, units, encoded_values) -> None:
+    # Regression test for GitHub issues #8271 and #9488
     var = Variable(
         ["time"],
         times,
@@ -1400,7 +1422,7 @@ def test_roundtrip_float_times() -> None:
     )
 
     encoded_var = conventions.encode_cf_variable(var)
-    np.testing.assert_array_equal(encoded_var, np.array([3653, 3653.25, 20.0]))
+    np.testing.assert_array_equal(encoded_var, encoded_values)
     assert encoded_var.attrs["units"] == units
     assert encoded_var.attrs["_FillValue"] == fill_value