Accept coordinates with MultiIndex (solve issue pydata#3008)

Author: ghislainp

xarray/core/coordinates.py

@@ -13,6 +13,7 @@
     cast,
 )
 
+import numpy as np
 import pandas as pd
 
 from . import formatting, indexing
@@ -106,9 +107,47 @@ def to_index(self, ordered_dims: Sequence[Hashable] = None) -> pd.Index:
             (dim,) = ordered_dims
             return self._data.get_index(dim)  # type: ignore
         else:
+            from pandas.core.arrays.categorical import factorize_from_iterable
+
             indexes = [self._data.get_index(k) for k in ordered_dims]  # type: ignore
-            names = list(ordered_dims)
-            return pd.MultiIndex.from_product(indexes, names=names)
+
+            # compute the sizes of the repeat and tile for the cartesian product
+            # (taken from pandas.core.reshape.util)
+            lenX = np.fromiter((len(index) for index in indexes), dtype=np.intp)
+            cumprodX = np.cumproduct(lenX)
+
+            if cumprodX[-1] != 0:
+                # sizes of the repeats
+                b = cumprodX[-1] / cumprodX
+            else:
+                # if any factor is empty, the cartesian product is empty
+                b = np.zeros_like(cumprodX)
+
+            # sizes of the tiles
+            a = np.roll(cumprodX, 1)
+            a[0] = 1
+
+            # loop over the indexes
+            # for each MultiIndex or Index compute the cartesian product of the codes
+
+            code_list = []
+            level_list = []
+            names = []
+
+            for i, index in enumerate(indexes):
+                if isinstance(index, pd.MultiIndex):
+                    codes, levels = index.codes, index.levels
+                else:
+                    code, level = factorize_from_iterable(index)
+                    codes = [code]
+                    levels = [level]
+                    
+                # compute the cartesian product
+                code_list += [np.tile(np.repeat(code, b[i]), a[i]) for code in codes]
+                level_list += levels
+                names += index.names
+
+            return pd.MultiIndex(level_list, code_list, names=names)
 
     def update(self, other: Mapping[Hashable, Any]) -> None:
         other_vars = getattr(other, "variables", other)