Repository: pydata/xarray

Issue:
test failure with upstream-dev
https://dev.azure.com/xarray/xarray/_build/results?buildId=1101&view=logs

```
=================================== FAILURES ===================================
_________________________ test_datetime_reduce[False] __________________________

dask = False

    @arm_xfail
    @pytest.mark.parametrize("dask", [False, True])
    def test_datetime_reduce(dask):
        time = np.array(pd.date_range("15/12/1999", periods=11))
        time[8:11] = np.nan
        da = DataArray(np.linspace(0, 365, num=11), dims="time", coords={"time": time})
    
        if dask and has_dask:
            chunks = {"time": 5}
            da = da.chunk(chunks)
    
        actual = da["time"].mean()
>       assert not pd.isnull(actual)
E       AssertionError: assert not True
E        +  where True = <function isna at 0x7f2475449ae8>(<xarray.DataArray 'time' ()>\narray('NaT', dtype='datetime64[ns]'))
E        +    where <function isna at 0x7f2475449ae8> = pd.isnull

xarray/tests/test_duck_array_ops.py:288: AssertionError
__________________________ test_datetime_reduce[True] __________________________

dask = True

    @arm_xfail
    @pytest.mark.parametrize("dask", [False, True])
    def test_datetime_reduce(dask):
        time = np.array(pd.date_range("15/12/1999", periods=11))
        time[8:11] = np.nan
        da = DataArray(np.linspace(0, 365, num=11), dims="time", coords={"time": time})
    
        if dask and has_dask:
            chunks = {"time": 5}
            da = da.chunk(chunks)
    
        actual = da["time"].mean()
>       assert not pd.isnull(actual)
E       AssertionError: assert not True
E        +  where True = <function isna at 0x7f2475449ae8>(<xarray.DataArray 'time' ()>\narray('NaT', dtype='datetime64[ns]'))
E        +    where <function isna at 0x7f2475449ae8> = pd.isnull

xarray/tests/test_duck_array_ops.py:288: AssertionError
=============================== warnings summary 
```

Discussion:
crusaderky: xref #3434
crusaderky: Ready for review and merge
mathause: Given https://github.com/numpy/numpy/pull/14841 - is this worth the extra complexity?

Currently xarray is only broken wrt numpy master and there is a chance https://github.com/numpy/numpy/pull/14841 gets into 1.18 - then we could just use `np.nanmin`?

Also - is this the only place of `min` for datetimes? Should this also go into the operator dispatcher?
dopplershift: IMO, it's always best to release code that you *know* will work rather than rely on upstream to get something into the next release in order for you not to be broken. I say that both as a downstream user of xarray and a library maintainers.
dcherian: LGTM. I think we should merge so we can get a release out.
crusaderky: @dopplershift fully agree.

> is this the only place of min for datetimes? 

This PR *does not* fix the nan-reductions for datetime arrays with NaT.

It *very narrowly* fixes the use case of **mean() and nanmean() with a numpy backend**, so that they work on numpy 1.18 as they did on 1.17. As a matter of fact they never worked with dask, were never tested, and still won't work after the PR.

To emphasize: I am not introducing any new feature or fixing bugs on numpy 1.17; I'm just getting rid of regressions on the upgrade from 1.17 to 1.18.

Specifically, the PR fixes this failing test:
https://github.com/pydata/xarray/blob/aa876cfd6b3b97ee5028d089ec741d057e3ed688/xarray/tests/test_duck_array_ops.py#L276-L293

which was not well formulated - despite its name, it doesn't test all reductions, only mean(), and in fact never runs on dask (the index of a dask array is not chunked).
I have now rewritten it  for the sake of clarity.

*If* numpy/numpy#14841 makes into master before 1.18 is published, and *if* it actually delivers on making reductions work with NaT (I did not test it yet), I'll be happy to write a second PR and rework everything.
Even then, a switch of this sort will still be needed:

```python
def datetime_nanmin(array, **kwargs):
    if LooseVersion(numpy.__version__) < "1.18":
        return np.min(array, **kwargs)
    else:
        return np.nanmin(array, **kwargs)
```
The actual code will be in fact more complicated because you also need to think about dask.

max-sixty: Agree—let's merge?
mathause: Ok, thanks for the clarification. 
crusaderky: Updated PR description:

Fix mean() and nanmean() for datetime64 arrays on numpy backend when upgrading from numpy 1.17 to 1.18.
All other nan-reductions on datetime64s were broken before and remain broken.
All reductions on datetime64s, including mean() and nanmean(), were broken on dask and remain broken.
crusaderky:  I just realised there are no tests at all for ``datetime_to_numeric(offset=None) `` on a dask backend. I suspect I broke something - I need to look at it better. 
pep8speaks: Hello @crusaderky! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`xarray/core/duck_array_ops.py`](https://github.com/pydata/xarray/blob/341fad3f84e5625629123fbb7cfd96b9c2d37d1d/xarray/core/duck_array_ops.py):

> [Line 13:1](https://github.com/pydata/xarray/blob/341fad3f84e5625629123fbb7cfd96b9c2d37d1d/xarray/core/duck_array_ops.py#L13): [F401](https://duckduckgo.com/?q=pep8%20F401) 'distutils.version.LooseVersion' imported but unused




crusaderky: ``datetime_to_numeric(offset=None)`` never worked with dask to begin with - there's a big fat TODO in it about it. I rewrote most of the PR; now mean() works with dask.

Pre-empting an observation: *yes*, the nanmin function I wrote *could* be used to work around the shortcomings of both numpy and dask on datetime64. I don't think xarray should do it and IMHO the workaround currently implemented for mean() should never have been done to begin with.

Again, this PR only fixes the regression on numpy 1.18.

crusaderky: Ready for review and merge

Patch:
diff --git a/ci/azure/install.yml b/ci/azure/install.yml
index fee886ba..baa69bcc 100644
--- a/ci/azure/install.yml
+++ b/ci/azure/install.yml
@@ -16,9 +16,9 @@ steps:
         --pre \
         --upgrade \
         matplotlib \
+        numpy \
         pandas \
         scipy
-        # numpy \  # FIXME https://github.com/pydata/xarray/issues/3409
     pip install \
         --no-deps \
         --upgrade \
diff --git a/ci/requirements/py36.yml b/ci/requirements/py36.yml
index 10fe6925..820160b1 100644
--- a/ci/requirements/py36.yml
+++ b/ci/requirements/py36.yml
@@ -25,7 +25,7 @@ dependencies:
   - nc-time-axis
   - netcdf4
   - numba
-  - numpy<1.18  # FIXME https://github.com/pydata/xarray/issues/3409
+  - numpy
   - pandas
   - pint
   - pip
diff --git a/ci/requirements/py37.yml b/ci/requirements/py37.yml
index 827c664a..4a7aaf7d 100644
--- a/ci/requirements/py37.yml
+++ b/ci/requirements/py37.yml
@@ -25,7 +25,7 @@ dependencies:
   - nc-time-axis
   - netcdf4
   - numba
-  - numpy<1.18  # FIXME https://github.com/pydata/xarray/issues/3409
+  - numpy
   - pandas
   - pint
   - pip
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index cb274bca..0c929b5b 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -115,6 +115,12 @@ Bug fixes
   (:issue:`3402`). By `Deepak Cherian <https://github.com/dcherian/>`_
 - Allow appending datetime and bool data variables to zarr stores.
   (:issue:`3480`). By `Akihiro Matsukawa <https://github.com/amatsukawa/>`_.
+- Add support for numpy >=1.18 (); bugfix mean() on datetime64 arrays on dask backend
+  (:issue:`3409`, :pull:`3537`). By `Guido Imperiale <https://github.com/crusaderky>`_.
+- Add support for pandas >=0.26 (:issue:`3440`).
+  By `Deepak Cherian <https://github.com/dcherian>`_.
+- Add support for pseudonetcdf >=3.1 (:pull:`3485`).
+  By `Barron Henderson <https://github.com/barronh>`_.
 
 Documentation
 ~~~~~~~~~~~~~
@@ -133,7 +139,6 @@ Documentation
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
-
 - Added integration tests against `pint <https://pint.readthedocs.io/>`_.
   (:pull:`3238`, :pull:`3447`, :pull:`3493`, :pull:`3508`)
   by `Justus Magin <https://github.com/keewis>`_.
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 5de25461..c631a4c1 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -5316,7 +5316,9 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
                 datetime_unit, _ = np.datetime_data(coord_var.dtype)
             elif datetime_unit is None:
                 datetime_unit = "s"  # Default to seconds for cftime objects
-            coord_var = datetime_to_numeric(coord_var, datetime_unit=datetime_unit)
+            coord_var = coord_var._replace(
+                data=datetime_to_numeric(coord_var.data, datetime_unit=datetime_unit)
+            )
 
         variables = {}
         coord_names = set()
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 71e79335..cf616acb 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -351,6 +351,26 @@ cumsum_1d.numeric_only = True
 _mean = _create_nan_agg_method("mean")
 
 
+def _datetime_nanmin(array):
+    """nanmin() function for datetime64.
+
+    Caveats that this function deals with:
+
+    - In numpy < 1.18, min() on datetime64 incorrectly ignores NaT
+    - numpy nanmin() don't work on datetime64 (all versions at the moment of writing)
+    - dask min() does not work on datetime64 (all versions at the moment of writing)
+    """
+    assert array.dtype.kind in "mM"
+    dtype = array.dtype
+    # (NaT).astype(float) does not produce NaN...
+    array = where(pandas_isnull(array), np.nan, array.astype(float))
+    array = min(array, skipna=True)
+    if isinstance(array, float):
+        array = np.array(array)
+    # ...but (NaN).astype("M8") does produce NaT
+    return array.astype(dtype)
+
+
 def datetime_to_numeric(array, offset=None, datetime_unit=None, dtype=float):
     """Convert an array containing datetime-like data to an array of floats.
 
@@ -370,7 +390,10 @@ def datetime_to_numeric(array, offset=None, datetime_unit=None, dtype=float):
     """
     # TODO: make this function dask-compatible?
     if offset is None:
-        offset = array.min()
+        if array.dtype.kind in "Mm":
+            offset = _datetime_nanmin(array)
+        else:
+            offset = min(array)
     array = array - offset
 
     if not hasattr(array, "dtype"):  # scalar is converted to 0d-array
@@ -401,7 +424,8 @@ def mean(array, axis=None, skipna=None, **kwargs):
 
     array = asarray(array)
     if array.dtype.kind in "Mm":
-        offset = min(array)
+        offset = _datetime_nanmin(array)
+
         # xarray always uses np.datetime64[ns] for np.datetime64 data
         dtype = "timedelta64[ns]"
         return (


Test patch:
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index be40ce7c..de074da5 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -5874,7 +5874,9 @@ def test_trapz_datetime(dask, which_datetime):
 
     actual = da.integrate("time", datetime_unit="D")
     expected_data = np.trapz(
-        da, duck_array_ops.datetime_to_numeric(da["time"], datetime_unit="D"), axis=0
+        da.data,
+        duck_array_ops.datetime_to_numeric(da["time"].data, datetime_unit="D"),
+        axis=0,
     )
     expected = xr.DataArray(
         expected_data,
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index f678af2f..aee7bbd6 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -274,23 +274,39 @@ def assert_dask_array(da, dask):
 
 
 @arm_xfail
-@pytest.mark.parametrize("dask", [False, True])
-def test_datetime_reduce(dask):
-    time = np.array(pd.date_range("15/12/1999", periods=11))
-    time[8:11] = np.nan
-    da = DataArray(np.linspace(0, 365, num=11), dims="time", coords={"time": time})
-
-    if dask and has_dask:
-        chunks = {"time": 5}
-        da = da.chunk(chunks)
-
-    actual = da["time"].mean()
-    assert not pd.isnull(actual)
-    actual = da["time"].mean(skipna=False)
-    assert pd.isnull(actual)
-
-    # test for a 0d array
-    assert da["time"][0].mean() == da["time"][:1].mean()
+@pytest.mark.parametrize("dask", [False, True] if has_dask else [False])
+def test_datetime_mean(dask):
+    # Note: only testing numpy, as dask is broken upstream
+    da = DataArray(
+        np.array(["2010-01-01", "NaT", "2010-01-03", "NaT", "NaT"], dtype="M8"),
+        dims=["time"],
+    )
+    if dask:
+        # Trigger use case where a chunk is full of NaT
+        da = da.chunk({"time": 3})
+
+    expect = DataArray(np.array("2010-01-02", dtype="M8"))
+    expect_nat = DataArray(np.array("NaT", dtype="M8"))
+
+    actual = da.mean()
+    if dask:
+        assert actual.chunks is not None
+    assert_equal(actual, expect)
+
+    actual = da.mean(skipna=False)
+    if dask:
+        assert actual.chunks is not None
+    assert_equal(actual, expect_nat)
+
+    # tests for 1d array full of NaT
+    assert_equal(da[[1]].mean(), expect_nat)
+    assert_equal(da[[1]].mean(skipna=False), expect_nat)
+
+    # tests for a 0d array
+    assert_equal(da[0].mean(), da[0])
+    assert_equal(da[0].mean(skipna=False), da[0])
+    assert_equal(da[1].mean(), expect_nat)
+    assert_equal(da[1].mean(skipna=False), expect_nat)
 
 
 @requires_cftime