Repository: pydata/xarray

Issue:
groupby_bins incorrect correspondence with labels
### What happened?

As of version 2023.4.0 groupby_bins stopped giving correct results. This was reported in #7751. The fix applied in 2023.4.1 changes the behaviour but does not fix the problem in my use case where strings are used as labels for the groups. It appears that the sorting that was implemented for the fix sorts based on name (so e.g. 'one' comes after 'four', see example) instead of on the lower edge. I could be wrong about the reason but accessing groups like data[group['label']] breaks 

### What did you expect to happen?

Until 2023.3.0 accessing groups by name gave the correct subset from the DataArray

### Minimal Complete Verifiable Example

```Python
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import sys

print(f"numpy version: {np.__version__}")
print(f"xarray version: {xr.__version__}")
print(f"pandas version: {pd.__version__}")
print(f"python version: {sys.version}")


# Generate random data 
# Make the coordiantes follow a normal distribution
np.random.seed(42)
coords = np.random.normal(5, 5, 1000)
bins = np.logspace(-4, 1, 10)
labels = ['one', 'two', 'three', 'four', 'five', 'six', 'seven',
          'eight', 'nine']
# xArray
# Make a mock dataarray
darr = xr.DataArray(coords, coords=[coords], dims=["coords"])
groups_xr = darr.groupby_bins("coords", bins, labels=labels).groups
print(groups_xr.keys())
for lab in labels:
    print(lab, end=' ')
    try:
        print(int(darr[groups_xr[lab]].count()))
    except:
        print("nan")
```


### MVCE confirmation

- [x] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
- [x] Complete example — the example is self-contained, including all data and the text of any traceback.
- [X] 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.
- [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

### Relevant log output

_No response_

### Anything else we need to know?

Output from 2023.3.0 (keys in random order, but correct access by name)

numpy version: 1.23.5
xarray version: 2023.3.0
pandas version: 1.5.3
python version: 3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]
dict_keys(['nine', 'eight', 'seven', 'five', 'six', 'four'])
one nan
two nan
three nan
four 1
five 2
six 9
seven 27
eight 153
nine 506

Output from 2023.4.1 - Keys alphabetically sorted, wrong access by name
numpy version: 1.23.5
xarray version: 2023.4.1
pandas version: 1.5.3
python version: 3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]
dict_keys(['eight', 'five', 'four', 'nine', 'seven', 'six'])
one nan
two nan
three nan
four 9
five 2
six 506
seven 153
eight 1
nine 27




### Environment

<details>

INSTALLED VERSIONS
------------------
commit: None
python: 3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]
python-bits: 64
OS: Linux
OS-release: 5.19.0-40-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: en_US.UTF-8
LANG: el_GR.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.10.7
libnetcdf: 4.8.1

xarray: 2023.4.1
pandas: 1.5.3
numpy: 1.23.5
scipy: 1.10.2.dev0+2411.f2ee108
netCDF4: 1.5.8
pydap: None
h5netcdf: 1.1.0
h5py: 3.6.0
Nio: None
zarr: 2.14.2
cftime: 1.5.2
nc_time_axis: None
PseudoNetCDF: None
iris: None
bottleneck: 1.3.7
dask: 2023.4.0
distributed: 2023.4.0
matplotlib: 3.7.1
cartopy: 0.21.1
seaborn: 0.12.2
numbagg: None
fsspec: 2023.4.0
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 67.6.1
pip: 23.1
conda: None
pytest: 6.2.5
mypy: 0.942
IPython: 8.12.0
sphinx: 4.3.2
</details>


Discussion:
gsieros: Just confirmed that the fix works for the actual application that I was using (where N-NW etc were used as labels in binning wind directions), so all good in that respect.
dcherian: Thanks @gsieros, and apologies for the trouble. Clearly our test suite was insufficient.

 I'll push out a release soon.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index cad27670..a393a114 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -35,6 +35,9 @@ Deprecations
 Bug fixes
 ~~~~~~~~~
 
+- Fix binning when ``labels`` is specified. (:issue:`7766`).
+  By `Deepak Cherian <https://github.com/dcherian>`_.
+
 
 Documentation
 ~~~~~~~~~~~~~
diff --git a/xarray/core/groupby.py b/xarray/core/groupby.py
index b786fa60..a6d63b7e 100644
--- a/xarray/core/groupby.py
+++ b/xarray/core/groupby.py
@@ -338,7 +338,8 @@ def _factorize_bins(
     if (codes == -1).all():
         raise ValueError(f"None of the data falls within bins with edges {bins!r}")
     full_index = binned.categories
-    unique_values = np.sort(binned.unique().dropna())
+    uniques = np.sort(pd.unique(codes))
+    unique_values = full_index[uniques[uniques != -1]]
     group_indices = [g for g in _codes_to_groups(codes, len(full_index)) if g]
 
     if len(group_indices) == 0:


Test patch:
diff --git a/xarray/tests/test_groupby.py b/xarray/tests/test_groupby.py
index 13e26954..4488bda3 100644
--- a/xarray/tests/test_groupby.py
+++ b/xarray/tests/test_groupby.py
@@ -1373,7 +1373,22 @@ class TestDataArrayGroupBy:
 
     @pytest.mark.parametrize("use_flox", [True, False])
     @pytest.mark.parametrize("coords", [np.arange(4), np.arange(4)[::-1], [2, 0, 3, 1]])
-    def test_groupby_bins(self, coords: np.typing.ArrayLike, use_flox: bool) -> None:
+    @pytest.mark.parametrize(
+        "cut_kwargs",
+        (
+            {"labels": None, "include_lowest": True},
+            {"labels": None, "include_lowest": False},
+            {"labels": ["a", "b"]},
+            {"labels": [1.2, 3.5]},
+            {"labels": ["b", "a"]},
+        ),
+    )
+    def test_groupby_bins(
+        self,
+        coords: np.typing.ArrayLike,
+        use_flox: bool,
+        cut_kwargs: dict,
+    ) -> None:
         array = DataArray(
             np.arange(4), dims="dim_0", coords={"dim_0": coords}, name="a"
         )
@@ -1384,11 +1399,10 @@ class TestDataArrayGroupBy:
         bins = [0, 1.5, 5]
 
         df = array.to_dataframe()
-        df["dim_0_bins"] = pd.cut(array["dim_0"], bins)
+        df["dim_0_bins"] = pd.cut(array["dim_0"], bins, **cut_kwargs)
 
         expected_df = df.groupby("dim_0_bins").sum()
         # TODO: can't convert df with IntervalIndex to Xarray
-
         expected = (
             expected_df.reset_index(drop=True)
             .to_xarray()
@@ -1397,25 +1411,55 @@ class TestDataArrayGroupBy:
         )
 
         with xr.set_options(use_flox=use_flox):
-            actual = array.groupby_bins("dim_0", bins=bins).sum()
+            actual = array.groupby_bins("dim_0", bins=bins, **cut_kwargs).sum()
             assert_identical(expected, actual)
 
-            actual = array.groupby_bins("dim_0", bins=bins, labels=[1.2, 3.5]).sum()
-            assert_identical(expected.assign_coords(dim_0_bins=[1.2, 3.5]), actual)
-
-            actual = array.groupby_bins("dim_0", bins=bins).map(lambda x: x.sum())
+            actual = array.groupby_bins("dim_0", bins=bins, **cut_kwargs).map(
+                lambda x: x.sum()
+            )
             assert_identical(expected, actual)
 
             # make sure original array dims are unchanged
             assert len(array.dim_0) == 4
 
-            da = xr.DataArray(np.ones((2, 3, 4)))
-            bins = [-1, 0, 1, 2]
-            with xr.set_options(use_flox=False):
-                actual = da.groupby_bins("dim_0", bins).mean(...)
-            with xr.set_options(use_flox=True):
-                expected = da.groupby_bins("dim_0", bins).mean(...)
-            assert_allclose(actual, expected)
+    def test_groupby_bins_ellipsis(self):
+        da = xr.DataArray(np.ones((2, 3, 4)))
+        bins = [-1, 0, 1, 2]
+        with xr.set_options(use_flox=False):
+            actual = da.groupby_bins("dim_0", bins).mean(...)
+        with xr.set_options(use_flox=True):
+            expected = da.groupby_bins("dim_0", bins).mean(...)
+        assert_allclose(actual, expected)
+
+    @pytest.mark.parametrize("use_flox", [True, False])
+    def test_groupby_bins_gives_correct_subset(self, use_flox: bool) -> None:
+        # GH7766
+        rng = np.random.default_rng(42)
+        coords = rng.normal(5, 5, 1000)
+        bins = np.logspace(-4, 1, 10)
+        labels = [
+            "one",
+            "two",
+            "three",
+            "four",
+            "five",
+            "six",
+            "seven",
+            "eight",
+            "nine",
+        ]
+        # xArray
+        # Make a mock dataarray
+        darr = xr.DataArray(coords, coords=[coords], dims=["coords"])
+        expected = xr.DataArray(
+            [np.nan, np.nan, 1, 1, 1, 8, 31, 104, 542],
+            dims="coords_bins",
+            coords={"coords_bins": labels},
+        )
+        gb = darr.groupby_bins("coords", bins, labels=labels)
+        with xr.set_options(use_flox=use_flox):
+            actual = gb.count()
+        assert_identical(actual, expected)
 
     def test_groupby_bins_empty(self):
         array = DataArray(np.arange(4), [("x", range(4))])