from __future__ import annotations import enum import functools import operator from collections import Counter, defaultdict from contextlib import suppress from dataclasses import dataclass, field from datetime import timedelta from html import escape from typing import TYPE_CHECKING, Any, Callable, Hashable, Iterable, Mapping import numpy as np import pandas as pd from packaging.version import Version from xarray.core import duck_array_ops from xarray.core.nputils import NumpyVIndexAdapter from xarray.core.options import OPTIONS from xarray.core.pycompat import ( array_type, integer_types, is_duck_dask_array, mod_version, ) from xarray.core.types import T_Xarray from xarray.core.utils import ( NDArrayMixin, either_dict_or_kwargs, get_valid_numpy_dtype, to_0d_array, ) if TYPE_CHECKING: from numpy.typing import DTypeLike from xarray.core.indexes import Index from xarray.core.variable import Variable @dataclass class IndexSelResult: """Index query results. Attributes ---------- dim_indexers: dict A dictionary where keys are array dimensions and values are location-based indexers. indexes: dict, optional New indexes to replace in the resulting DataArray or Dataset. variables : dict, optional New variables to replace in the resulting DataArray or Dataset. drop_coords : list, optional Coordinate(s) to drop in the resulting DataArray or Dataset. drop_indexes : list, optional Index(es) to drop in the resulting DataArray or Dataset. rename_dims : dict, optional A dictionary in the form ``{old_dim: new_dim}`` for dimension(s) to rename in the resulting DataArray or Dataset. """ dim_indexers: dict[Any, Any] indexes: dict[Any, Index] = field(default_factory=dict) variables: dict[Any, Variable] = field(default_factory=dict) drop_coords: list[Hashable] = field(default_factory=list) drop_indexes: list[Hashable] = field(default_factory=list) rename_dims: dict[Any, Hashable] = field(default_factory=dict) def as_tuple(self): """Unlike ``dataclasses.astuple``, return a shallow copy. See https://stackoverflow.com/a/51802661 """ return ( self.dim_indexers, self.indexes, self.variables, self.drop_coords, self.drop_indexes, self.rename_dims, ) def merge_sel_results(results: list[IndexSelResult]) -> IndexSelResult: all_dims_count = Counter([dim for res in results for dim in res.dim_indexers]) duplicate_dims = {k: v for k, v in all_dims_count.items() if v > 1} if duplicate_dims: # TODO: this message is not right when combining indexe(s) queries with # location-based indexing on a dimension with no dimension-coordinate (failback) fmt_dims = [ f"{dim!r}: {count} indexes involved" for dim, count in duplicate_dims.items() ] raise ValueError( "Xarray does not support label-based selection with more than one index " "over the following dimension(s):\n" + "\n".join(fmt_dims) + "\nSuggestion: use a multi-index for each of those dimension(s)." ) dim_indexers = {} indexes = {} variables = {} drop_coords = [] drop_indexes = [] rename_dims = {} for res in results: dim_indexers.update(res.dim_indexers) indexes.update(res.indexes) variables.update(res.variables) drop_coords += res.drop_coords drop_indexes += res.drop_indexes rename_dims.update(res.rename_dims) return IndexSelResult( dim_indexers, indexes, variables, drop_coords, drop_indexes, rename_dims ) def group_indexers_by_index( obj: T_Xarray, indexers: Mapping[Any, Any], options: Mapping[str, Any], ) -> list[tuple[Index, dict[Any, Any]]]: """Returns a list of unique indexes and their corresponding indexers.""" unique_indexes = {} grouped_indexers: Mapping[int | None, dict] = defaultdict(dict) for key, label in indexers.items(): index: Index = obj.xindexes.get(key, None) if index is not None: index_id = id(index) unique_indexes[index_id] = index grouped_indexers[index_id][key] = label elif key in obj.coords: raise KeyError(f"no index found for coordinate {key!r}") elif key not in obj.dims: raise KeyError(f"{key!r} is not a valid dimension or coordinate") elif len(options): raise ValueError( f"cannot supply selection options {options!r} for dimension {key!r}" "that has no associated coordinate or index" ) else: # key is a dimension without a "dimension-coordinate" # failback to location-based selection # TODO: depreciate this implicit behavior and suggest using isel instead? unique_indexes[None] = None grouped_indexers[None][key] = label return [(unique_indexes[k], grouped_indexers[k]) for k in unique_indexes] def map_index_queries( obj: T_Xarray, indexers: Mapping[Any, Any], method=None, tolerance=None, **indexers_kwargs: Any, ) -> IndexSelResult: """Execute index queries from a DataArray / Dataset and label-based indexers and return the (merged) query results. """ from xarray.core.dataarray import DataArray # TODO benbovy - flexible indexes: remove when custom index options are available if method is None and tolerance is None: options = {} else: options = {"method": method, "tolerance": tolerance} indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "map_index_queries") grouped_indexers = group_indexers_by_index(obj, indexers, options) results = [] for index, labels in grouped_indexers: if index is None: # forward dimension indexers with no index/coordinate results.append(IndexSelResult(labels)) else: results.append(index.sel(labels, **options)) merged = merge_sel_results(results) # drop dimension coordinates found in dimension indexers # (also drop multi-index if any) # (.sel() already ensures alignment) for k, v in merged.dim_indexers.items(): if isinstance(v, DataArray): if k in v._indexes: v = v.reset_index(k) drop_coords = [name for name in v._coords if name in merged.dim_indexers] merged.dim_indexers[k] = v.drop_vars(drop_coords) return merged def expanded_indexer(key, ndim): """Given a key for indexing an ndarray, return an equivalent key which is a tuple with length equal to the number of dimensions. The expansion is done by replacing all `Ellipsis` items with the right number of full slices and then padding the key with full slices so that it reaches the appropriate dimensionality. """ if not isinstance(key, tuple): # numpy treats non-tuple keys equivalent to tuples of length 1 key = (key,) new_key = [] # handling Ellipsis right is a little tricky, see: # https://numpy.org/doc/stable/reference/arrays.indexing.html#advanced-indexing found_ellipsis = False for k in key: if k is Ellipsis: if not found_ellipsis: new_key.extend((ndim + 1 - len(key)) * [slice(None)]) found_ellipsis = True else: new_key.append(slice(None)) else: new_key.append(k) if len(new_key) > ndim: raise IndexError("too many indices") new_key.extend((ndim - len(new_key)) * [slice(None)]) return tuple(new_key) def _expand_slice(slice_, size): return np.arange(*slice_.indices(size)) def _normalize_slice(sl, size): """Ensure that given slice only contains positive start and stop values (stop can be -1 for full-size slices with negative steps, e.g. [-10::-1])""" return slice(*sl.indices(size)) def slice_slice(old_slice, applied_slice, size): """Given a slice and the size of the dimension to which it will be applied, index it with another slice to return a new slice equivalent to applying the slices sequentially """ old_slice = _normalize_slice(old_slice, size) size_after_old_slice = len(range(old_slice.start, old_slice.stop, old_slice.step)) if size_after_old_slice == 0: # nothing left after applying first slice return slice(0) applied_slice = _normalize_slice(applied_slice, size_after_old_slice) start = old_slice.start + applied_slice.start * old_slice.step if start < 0: # nothing left after applying second slice # (can only happen for old_slice.step < 0, e.g. [10::-1], [20:]) return slice(0) stop = old_slice.start + applied_slice.stop * old_slice.step if stop < 0: stop = None step = old_slice.step * applied_slice.step return slice(start, stop, step) def _index_indexer_1d(old_indexer, applied_indexer, size): assert isinstance(applied_indexer, integer_types + (slice, np.ndarray)) if isinstance(applied_indexer, slice) and applied_indexer == slice(None): # shortcut for the usual case return old_indexer if isinstance(old_indexer, slice): if isinstance(applied_indexer, slice): indexer = slice_slice(old_indexer, applied_indexer, size) else: indexer = _expand_slice(old_indexer, size)[applied_indexer] else: indexer = old_indexer[applied_indexer] return indexer class ExplicitIndexer: """Base class for explicit indexer objects. ExplicitIndexer objects wrap a tuple of values given by their ``tuple`` property. These tuples should always have length equal to the number of dimensions on the indexed array. Do not instantiate BaseIndexer objects directly: instead, use one of the sub-classes BasicIndexer, OuterIndexer or VectorizedIndexer. """ __slots__ = ("_key",) def __init__(self, key): if type(self) is ExplicitIndexer: raise TypeError("cannot instantiate base ExplicitIndexer objects") self._key = tuple(key) @property def tuple(self): return self._key def __repr__(self): return f"{type(self).__name__}({self.tuple})" def as_integer_or_none(value): return None if value is None else operator.index(value) def as_integer_slice(value): start = as_integer_or_none(value.start) stop = as_integer_or_none(value.stop) step = as_integer_or_none(value.step) return slice(start, stop, step) class BasicIndexer(ExplicitIndexer): """Tuple for basic indexing. All elements should be int or slice objects. Indexing follows NumPy's rules for basic indexing: each axis is independently sliced and axes indexed with an integer are dropped from the result. """ __slots__ = () def __init__(self, key): if not isinstance(key, tuple): raise TypeError(f"key must be a tuple: {key!r}") new_key = [] for k in key: if isinstance(k, integer_types): k = int(k) elif isinstance(k, slice): k = as_integer_slice(k) else: raise TypeError( f"unexpected indexer type for {type(self).__name__}: {k!r}" ) new_key.append(k) super().__init__(new_key) class OuterIndexer(ExplicitIndexer): """Tuple for outer/orthogonal indexing. All elements should be int, slice or 1-dimensional np.ndarray objects with an integer dtype. Indexing is applied independently along each axis, and axes indexed with an integer are dropped from the result. This type of indexing works like MATLAB/Fortran. """ __slots__ = () def __init__(self, key): if not isinstance(key, tuple): raise TypeError(f"key must be a tuple: {key!r}") new_key = [] for k in key: if isinstance(k, integer_types): k = int(k) elif isinstance(k, slice): k = as_integer_slice(k) elif isinstance(k, np.ndarray): if not np.issubdtype(k.dtype, np.integer): raise TypeError( f"invalid indexer array, does not have integer dtype: {k!r}" ) if k.ndim != 1: raise TypeError( f"invalid indexer array for {type(self).__name__}; must have " f"exactly 1 dimension: {k!r}" ) k = np.asarray(k, dtype=np.int64) else: raise TypeError( f"unexpected indexer type for {type(self).__name__}: {k!r}" ) new_key.append(k) super().__init__(new_key) class VectorizedIndexer(ExplicitIndexer): """Tuple for vectorized indexing. All elements should be slice or N-dimensional np.ndarray objects with an integer dtype and the same number of dimensions. Indexing follows proposed rules for np.ndarray.vindex, which matches NumPy's advanced indexing rules (including broadcasting) except sliced axes are always moved to the end: https://github.com/numpy/numpy/pull/6256 """ __slots__ = () def __init__(self, key): if not isinstance(key, tuple): raise TypeError(f"key must be a tuple: {key!r}") new_key = [] ndim = None for k in key: if isinstance(k, slice): k = as_integer_slice(k) elif isinstance(k, np.ndarray): if not np.issubdtype(k.dtype, np.integer): raise TypeError( f"invalid indexer array, does not have integer dtype: {k!r}" ) if ndim is None: ndim = k.ndim elif ndim != k.ndim: ndims = [k.ndim for k in key if isinstance(k, np.ndarray)] raise ValueError( "invalid indexer key: ndarray arguments " f"have different numbers of dimensions: {ndims}" ) k = np.asarray(k, dtype=np.int64) else: raise TypeError( f"unexpected indexer type for {type(self).__name__}: {k!r}" ) new_key.append(k) super().__init__(new_key) class ExplicitlyIndexed: """Mixin to mark support for Indexer subclasses in indexing.""" __slots__ = () class ExplicitlyIndexedNDArrayMixin(NDArrayMixin, ExplicitlyIndexed): __slots__ = () def __array__(self, dtype=None): key = BasicIndexer((slice(None),) * self.ndim) return np.asarray(self[key], dtype=dtype) class ImplicitToExplicitIndexingAdapter(NDArrayMixin): """Wrap an array, converting tuples into the indicated explicit indexer.""" __slots__ = ("array", "indexer_cls") def __init__(self, array, indexer_cls=BasicIndexer): self.array = as_indexable(array) self.indexer_cls = indexer_cls def __array__(self, dtype=None): return np.asarray(self.array, dtype=dtype) def __getitem__(self, key): key = expanded_indexer(key, self.ndim) result = self.array[self.indexer_cls(key)] if isinstance(result, ExplicitlyIndexed): return type(self)(result, self.indexer_cls) else: # Sometimes explicitly indexed arrays return NumPy arrays or # scalars. return result class LazilyIndexedArray(ExplicitlyIndexedNDArrayMixin): """Wrap an array to make basic and outer indexing lazy.""" __slots__ = ("array", "key") def __init__(self, array, key=None): """ Parameters ---------- array : array_like Array like object to index. key : ExplicitIndexer, optional Array indexer. If provided, it is assumed to already be in canonical expanded form. """ if isinstance(array, type(self)) and key is None: # unwrap key = array.key array = array.array if key is None: key = BasicIndexer((slice(None),) * array.ndim) self.array = as_indexable(array) self.key = key def _updated_key(self, new_key): iter_new_key = iter(expanded_indexer(new_key.tuple, self.ndim)) full_key = [] for size, k in zip(self.array.shape, self.key.tuple): if isinstance(k, integer_types): full_key.append(k) else: full_key.append(_index_indexer_1d(k, next(iter_new_key), size)) full_key = tuple(full_key) if all(isinstance(k, integer_types + (slice,)) for k in full_key): return BasicIndexer(full_key) return OuterIndexer(full_key) @property def shape(self) -> tuple[int, ...]: shape = [] for size, k in zip(self.array.shape, self.key.tuple): if isinstance(k, slice): shape.append(len(range(*k.indices(size)))) elif isinstance(k, np.ndarray): shape.append(k.size) return tuple(shape) def __array__(self, dtype=None): array = as_indexable(self.array) return np.asarray(array[self.key], dtype=None) def transpose(self, order): return LazilyVectorizedIndexedArray(self.array, self.key).transpose(order) def __getitem__(self, indexer): if isinstance(indexer, VectorizedIndexer): array = LazilyVectorizedIndexedArray(self.array, self.key) return array[indexer] return type(self)(self.array, self._updated_key(indexer)) def __setitem__(self, key, value): if isinstance(key, VectorizedIndexer): raise NotImplementedError( "Lazy item assignment with the vectorized indexer is not yet " "implemented. Load your data first by .load() or compute()." ) full_key = self._updated_key(key) self.array[full_key] = value def __repr__(self): return f"{type(self).__name__}(array={self.array!r}, key={self.key!r})" # keep an alias to the old name for external backends pydata/xarray#5111 LazilyOuterIndexedArray = LazilyIndexedArray class LazilyVectorizedIndexedArray(ExplicitlyIndexedNDArrayMixin): """Wrap an array to make vectorized indexing lazy.""" __slots__ = ("array", "key") def __init__(self, array, key): """ Parameters ---------- array : array_like Array like object to index. key : VectorizedIndexer """ if isinstance(key, (BasicIndexer, OuterIndexer)): self.key = _outer_to_vectorized_indexer(key, array.shape) else: self.key = _arrayize_vectorized_indexer(key, array.shape) self.array = as_indexable(array) @property def shape(self) -> tuple[int, ...]: return np.broadcast(*self.key.tuple).shape def __array__(self, dtype=None): return np.asarray(self.array[self.key], dtype=None) def _updated_key(self, new_key): return _combine_indexers(self.key, self.shape, new_key) def __getitem__(self, indexer): # If the indexed array becomes a scalar, return LazilyIndexedArray if all(isinstance(ind, integer_types) for ind in indexer.tuple): key = BasicIndexer(tuple(k[indexer.tuple] for k in self.key.tuple)) return LazilyIndexedArray(self.array, key) return type(self)(self.array, self._updated_key(indexer)) def transpose(self, order): key = VectorizedIndexer(tuple(k.transpose(order) for k in self.key.tuple)) return type(self)(self.array, key) def __setitem__(self, key, value): raise NotImplementedError( "Lazy item assignment with the vectorized indexer is not yet " "implemented. Load your data first by .load() or compute()." ) def __repr__(self): return f"{type(self).__name__}(array={self.array!r}, key={self.key!r})" def _wrap_numpy_scalars(array): """Wrap NumPy scalars in 0d arrays.""" if np.isscalar(array): return np.array(array) else: return array class CopyOnWriteArray(ExplicitlyIndexedNDArrayMixin): __slots__ = ("array", "_copied") def __init__(self, array): self.array = as_indexable(array) self._copied = False def _ensure_copied(self): if not self._copied: self.array = as_indexable(np.array(self.array)) self._copied = True def __array__(self, dtype=None): return np.asarray(self.array, dtype=dtype) def __getitem__(self, key): return type(self)(_wrap_numpy_scalars(self.array[key])) def transpose(self, order): return self.array.transpose(order) def __setitem__(self, key, value): self._ensure_copied() self.array[key] = value def __deepcopy__(self, memo): # CopyOnWriteArray is used to wrap backend array objects, which might # point to files on disk, so we can't rely on the default deepcopy # implementation. return type(self)(self.array) class MemoryCachedArray(ExplicitlyIndexedNDArrayMixin): __slots__ = ("array",) def __init__(self, array): self.array = _wrap_numpy_scalars(as_indexable(array)) def _ensure_cached(self): if not isinstance(self.array, NumpyIndexingAdapter): self.array = NumpyIndexingAdapter(np.asarray(self.array)) def __array__(self, dtype=None): self._ensure_cached() return np.asarray(self.array, dtype=dtype) def __getitem__(self, key): return type(self)(_wrap_numpy_scalars(self.array[key])) def transpose(self, order): return self.array.transpose(order) def __setitem__(self, key, value): self.array[key] = value def as_indexable(array): """ This function always returns a ExplicitlyIndexed subclass, so that the vectorized indexing is always possible with the returned object. """ if isinstance(array, ExplicitlyIndexed): return array if isinstance(array, np.ndarray): return NumpyIndexingAdapter(array) if isinstance(array, pd.Index): return PandasIndexingAdapter(array) if is_duck_dask_array(array): return DaskIndexingAdapter(array) if hasattr(array, "__array_function__"): return NdArrayLikeIndexingAdapter(array) if hasattr(array, "__array_namespace__"): return ArrayApiIndexingAdapter(array) raise TypeError(f"Invalid array type: {type(array)}") def _outer_to_vectorized_indexer(key, shape): """Convert an OuterIndexer into an vectorized indexer. Parameters ---------- key : Outer/Basic Indexer An indexer to convert. shape : tuple Shape of the array subject to the indexing. Returns ------- VectorizedIndexer Tuple suitable for use to index a NumPy array with vectorized indexing. Each element is an array: broadcasting them together gives the shape of the result. """ key = key.tuple n_dim = len([k for k in key if not isinstance(k, integer_types)]) i_dim = 0 new_key = [] for k, size in zip(key, shape): if isinstance(k, integer_types): new_key.append(np.array(k).reshape((1,) * n_dim)) else: # np.ndarray or slice if isinstance(k, slice): k = np.arange(*k.indices(size)) assert k.dtype.kind in {"i", "u"} shape = [(1,) * i_dim + (k.size,) + (1,) * (n_dim - i_dim - 1)] new_key.append(k.reshape(*shape)) i_dim += 1 return VectorizedIndexer(tuple(new_key)) def _outer_to_numpy_indexer(key, shape): """Convert an OuterIndexer into an indexer for NumPy. Parameters ---------- key : Basic/OuterIndexer An indexer to convert. shape : tuple Shape of the array subject to the indexing. Returns ------- tuple Tuple suitable for use to index a NumPy array. """ if len([k for k in key.tuple if not isinstance(k, slice)]) <= 1: # If there is only one vector and all others are slice, # it can be safely used in mixed basic/advanced indexing. # Boolean index should already be converted to integer array. return key.tuple else: return _outer_to_vectorized_indexer(key, shape).tuple def _combine_indexers(old_key, shape, new_key): """Combine two indexers. Parameters ---------- old_key : ExplicitIndexer The first indexer for the original array shape : tuple of ints Shape of the original array to be indexed by old_key new_key The second indexer for indexing original[old_key] """ if not isinstance(old_key, VectorizedIndexer): old_key = _outer_to_vectorized_indexer(old_key, shape) if len(old_key.tuple) == 0: return new_key new_shape = np.broadcast(*old_key.tuple).shape if isinstance(new_key, VectorizedIndexer): new_key = _arrayize_vectorized_indexer(new_key, new_shape) else: new_key = _outer_to_vectorized_indexer(new_key, new_shape) return VectorizedIndexer( tuple(o[new_key.tuple] for o in np.broadcast_arrays(*old_key.tuple)) ) @enum.unique class IndexingSupport(enum.Enum): # for backends that support only basic indexer BASIC = 0 # for backends that support basic / outer indexer OUTER = 1 # for backends that support outer indexer including at most 1 vector. OUTER_1VECTOR = 2 # for backends that support full vectorized indexer. VECTORIZED = 3 def explicit_indexing_adapter( key: ExplicitIndexer, shape: tuple[int, ...], indexing_support: IndexingSupport, raw_indexing_method: Callable, ) -> Any: """Support explicit indexing by delegating to a raw indexing method. Outer and/or vectorized indexers are supported by indexing a second time with a NumPy array. Parameters ---------- key : ExplicitIndexer Explicit indexing object. shape : Tuple[int, ...] Shape of the indexed array. indexing_support : IndexingSupport enum Form of indexing supported by raw_indexing_method. raw_indexing_method : callable Function (like ndarray.__getitem__) that when called with indexing key in the form of a tuple returns an indexed array. Returns ------- Indexing result, in the form of a duck numpy-array. """ raw_key, numpy_indices = decompose_indexer(key, shape, indexing_support) result = raw_indexing_method(raw_key.tuple) if numpy_indices.tuple: # index the loaded np.ndarray result = NumpyIndexingAdapter(np.asarray(result))[numpy_indices] return result def decompose_indexer( indexer: ExplicitIndexer, shape: tuple[int, ...], indexing_support: IndexingSupport ) -> tuple[ExplicitIndexer, ExplicitIndexer]: if isinstance(indexer, VectorizedIndexer): return _decompose_vectorized_indexer(indexer, shape, indexing_support) if isinstance(indexer, (BasicIndexer, OuterIndexer)): return _decompose_outer_indexer(indexer, shape, indexing_support) raise TypeError(f"unexpected key type: {indexer}") def _decompose_slice(key, size): """convert a slice to successive two slices. The first slice always has a positive step. """ start, stop, step = key.indices(size) if step > 0: # If key already has a positive step, use it as is in the backend return key, slice(None) else: # determine stop precisely for step > 1 case # e.g. [98:2:-2] -> [98:3:-2] stop = start + int((stop - start - 1) / step) * step + 1 start, stop = stop + 1, start + 1 return slice(start, stop, -step), slice(None, None, -1) def _decompose_vectorized_indexer( indexer: VectorizedIndexer, shape: tuple[int, ...], indexing_support: IndexingSupport, ) -> tuple[ExplicitIndexer, ExplicitIndexer]: """ Decompose vectorized indexer to the successive two indexers, where the first indexer will be used to index backend arrays, while the second one is used to index loaded on-memory np.ndarray. Parameters ---------- indexer : VectorizedIndexer indexing_support : one of IndexerSupport entries Returns ------- backend_indexer: OuterIndexer or BasicIndexer np_indexers: an ExplicitIndexer (VectorizedIndexer / BasicIndexer) Notes ----- This function is used to realize the vectorized indexing for the backend arrays that only support basic or outer indexing. As an example, let us consider to index a few elements from a backend array with a vectorized indexer ([0, 3, 1], [2, 3, 2]). Even if the backend array only supports outer indexing, it is more efficient to load a subslice of the array than loading the entire array, >>> array = np.arange(36).reshape(6, 6) >>> backend_indexer = OuterIndexer((np.array([0, 1, 3]), np.array([2, 3]))) >>> # load subslice of the array ... array = NumpyIndexingAdapter(array)[backend_indexer] >>> np_indexer = VectorizedIndexer((np.array([0, 2, 1]), np.array([0, 1, 0]))) >>> # vectorized indexing for on-memory np.ndarray. ... NumpyIndexingAdapter(array)[np_indexer] array([ 2, 21, 8]) """ assert isinstance(indexer, VectorizedIndexer) if indexing_support is IndexingSupport.VECTORIZED: return indexer, BasicIndexer(()) backend_indexer_elems = [] np_indexer_elems = [] # convert negative indices indexer_elems = [ np.where(k < 0, k + s, k) if isinstance(k, np.ndarray) else k for k, s in zip(indexer.tuple, shape) ] for k, s in zip(indexer_elems, shape): if isinstance(k, slice): # If it is a slice, then we will slice it as-is # (but make its step positive) in the backend, # and then use all of it (slice(None)) for the in-memory portion. bk_slice, np_slice = _decompose_slice(k, s) backend_indexer_elems.append(bk_slice) np_indexer_elems.append(np_slice) else: # If it is a (multidimensional) np.ndarray, just pickup the used # keys without duplication and store them as a 1d-np.ndarray. oind, vind = np.unique(k, return_inverse=True) backend_indexer_elems.append(oind) np_indexer_elems.append(vind.reshape(*k.shape)) backend_indexer = OuterIndexer(tuple(backend_indexer_elems)) np_indexer = VectorizedIndexer(tuple(np_indexer_elems)) if indexing_support is IndexingSupport.OUTER: return backend_indexer, np_indexer # If the backend does not support outer indexing, # backend_indexer (OuterIndexer) is also decomposed. backend_indexer1, np_indexer1 = _decompose_outer_indexer( backend_indexer, shape, indexing_support ) np_indexer = _combine_indexers(np_indexer1, shape, np_indexer) return backend_indexer1, np_indexer def _decompose_outer_indexer( indexer: BasicIndexer | OuterIndexer, shape: tuple[int, ...], indexing_support: IndexingSupport, ) -> tuple[ExplicitIndexer, ExplicitIndexer]: """ Decompose outer indexer to the successive two indexers, where the first indexer will be used to index backend arrays, while the second one is used to index the loaded on-memory np.ndarray. Parameters ---------- indexer : OuterIndexer or BasicIndexer indexing_support : One of the entries of IndexingSupport Returns ------- backend_indexer: OuterIndexer or BasicIndexer np_indexers: an ExplicitIndexer (OuterIndexer / BasicIndexer) Notes ----- This function is used to realize the vectorized indexing for the backend arrays that only support basic or outer indexing. As an example, let us consider to index a few elements from a backend array with a orthogonal indexer ([0, 3, 1], [2, 3, 2]). Even if the backend array only supports basic indexing, it is more efficient to load a subslice of the array than loading the entire array, >>> array = np.arange(36).reshape(6, 6) >>> backend_indexer = BasicIndexer((slice(0, 3), slice(2, 4))) >>> # load subslice of the array ... array = NumpyIndexingAdapter(array)[backend_indexer] >>> np_indexer = OuterIndexer((np.array([0, 2, 1]), np.array([0, 1, 0]))) >>> # outer indexing for on-memory np.ndarray. ... NumpyIndexingAdapter(array)[np_indexer] array([[ 2, 3, 2], [14, 15, 14], [ 8, 9, 8]]) """ if indexing_support == IndexingSupport.VECTORIZED: return indexer, BasicIndexer(()) assert isinstance(indexer, (OuterIndexer, BasicIndexer)) backend_indexer: list[Any] = [] np_indexer = [] # make indexer positive pos_indexer: list[np.ndarray | int | np.number] = [] for k, s in zip(indexer.tuple, shape): if isinstance(k, np.ndarray): pos_indexer.append(np.where(k < 0, k + s, k)) elif isinstance(k, integer_types) and k < 0: pos_indexer.append(k + s) else: pos_indexer.append(k) indexer_elems = pos_indexer if indexing_support is IndexingSupport.OUTER_1VECTOR: # some backends such as h5py supports only 1 vector in indexers # We choose the most efficient axis gains = [ (np.max(k) - np.min(k) + 1.0) / len(np.unique(k)) if isinstance(k, np.ndarray) else 0 for k in indexer_elems ] array_index = np.argmax(np.array(gains)) if len(gains) > 0 else None for i, (k, s) in enumerate(zip(indexer_elems, shape)): if isinstance(k, np.ndarray) and i != array_index: # np.ndarray key is converted to slice that covers the entire # entries of this key. backend_indexer.append(slice(np.min(k), np.max(k) + 1)) np_indexer.append(k - np.min(k)) elif isinstance(k, np.ndarray): # Remove duplicates and sort them in the increasing order pkey, ekey = np.unique(k, return_inverse=True) backend_indexer.append(pkey) np_indexer.append(ekey) elif isinstance(k, integer_types): backend_indexer.append(k) else: # slice: convert positive step slice for backend bk_slice, np_slice = _decompose_slice(k, s) backend_indexer.append(bk_slice) np_indexer.append(np_slice) return (OuterIndexer(tuple(backend_indexer)), OuterIndexer(tuple(np_indexer))) if indexing_support == IndexingSupport.OUTER: for k, s in zip(indexer_elems, shape): if isinstance(k, slice): # slice: convert positive step slice for backend bk_slice, np_slice = _decompose_slice(k, s) backend_indexer.append(bk_slice) np_indexer.append(np_slice) elif isinstance(k, integer_types): backend_indexer.append(k) elif isinstance(k, np.ndarray) and (np.diff(k) >= 0).all(): backend_indexer.append(k) np_indexer.append(slice(None)) else: # Remove duplicates and sort them in the increasing order oind, vind = np.unique(k, return_inverse=True) backend_indexer.append(oind) np_indexer.append(vind.reshape(*k.shape)) return (OuterIndexer(tuple(backend_indexer)), OuterIndexer(tuple(np_indexer))) # basic indexer assert indexing_support == IndexingSupport.BASIC for k, s in zip(indexer_elems, shape): if isinstance(k, np.ndarray): # np.ndarray key is converted to slice that covers the entire # entries of this key. backend_indexer.append(slice(np.min(k), np.max(k) + 1)) np_indexer.append(k - np.min(k)) elif isinstance(k, integer_types): backend_indexer.append(k) else: # slice: convert positive step slice for backend bk_slice, np_slice = _decompose_slice(k, s) backend_indexer.append(bk_slice) np_indexer.append(np_slice) return (BasicIndexer(tuple(backend_indexer)), OuterIndexer(tuple(np_indexer))) def _arrayize_vectorized_indexer(indexer, shape): """Return an identical vindex but slices are replaced by arrays""" slices = [v for v in indexer.tuple if isinstance(v, slice)] if len(slices) == 0: return indexer arrays = [v for v in indexer.tuple if isinstance(v, np.ndarray)] n_dim = arrays[0].ndim if len(arrays) > 0 else 0 i_dim = 0 new_key = [] for v, size in zip(indexer.tuple, shape): if isinstance(v, np.ndarray): new_key.append(np.reshape(v, v.shape + (1,) * len(slices))) else: # slice shape = (1,) * (n_dim + i_dim) + (-1,) + (1,) * (len(slices) - i_dim - 1) new_key.append(np.arange(*v.indices(size)).reshape(shape)) i_dim += 1 return VectorizedIndexer(tuple(new_key)) def _dask_array_with_chunks_hint(array, chunks): """Create a dask array using the chunks hint for dimensions of size > 1.""" import dask.array as da if len(chunks) < array.ndim: raise ValueError("not enough chunks in hint") new_chunks = [] for chunk, size in zip(chunks, array.shape): new_chunks.append(chunk if size > 1 else (1,)) return da.from_array(array, new_chunks) def _logical_any(args): return functools.reduce(operator.or_, args) def _masked_result_drop_slice(key, data=None): key = (k for k in key if not isinstance(k, slice)) chunks_hint = getattr(data, "chunks", None) new_keys = [] for k in key: if isinstance(k, np.ndarray): if is_duck_dask_array(data): new_keys.append(_dask_array_with_chunks_hint(k, chunks_hint)) elif isinstance(data, array_type("sparse")): import sparse new_keys.append(sparse.COO.from_numpy(k)) else: new_keys.append(k) else: new_keys.append(k) mask = _logical_any(k == -1 for k in new_keys) return mask def create_mask(indexer, shape, data=None): """Create a mask for indexing with a fill-value. Parameters ---------- indexer : ExplicitIndexer Indexer with -1 in integer or ndarray value to indicate locations in the result that should be masked. shape : tuple Shape of the array being indexed. data : optional Data for which mask is being created. If data is a dask arrays, its chunks are used as a hint for chunks on the resulting mask. If data is a sparse array, the returned mask is also a sparse array. Returns ------- mask : bool, np.ndarray, SparseArray or dask.array.Array with dtype=bool Same type as data. Has the same shape as the indexing result. """ if isinstance(indexer, OuterIndexer): key = _outer_to_vectorized_indexer(indexer, shape).tuple assert not any(isinstance(k, slice) for k in key) mask = _masked_result_drop_slice(key, data) elif isinstance(indexer, VectorizedIndexer): key = indexer.tuple base_mask = _masked_result_drop_slice(key, data) slice_shape = tuple( np.arange(*k.indices(size)).size for k, size in zip(key, shape) if isinstance(k, slice) ) expanded_mask = base_mask[(Ellipsis,) + (np.newaxis,) * len(slice_shape)] mask = duck_array_ops.broadcast_to(expanded_mask, base_mask.shape + slice_shape) elif isinstance(indexer, BasicIndexer): mask = any(k == -1 for k in indexer.tuple) else: raise TypeError(f"unexpected key type: {type(indexer)}") return mask def _posify_mask_subindexer(index): """Convert masked indices in a flat array to the nearest unmasked index. Parameters ---------- index : np.ndarray One dimensional ndarray with dtype=int. Returns ------- np.ndarray One dimensional ndarray with all values equal to -1 replaced by an adjacent non-masked element. """ masked = index == -1 unmasked_locs = np.flatnonzero(~masked) if not unmasked_locs.size: # indexing unmasked_locs is invalid return np.zeros_like(index) masked_locs = np.flatnonzero(masked) prev_value = np.maximum(0, np.searchsorted(unmasked_locs, masked_locs) - 1) new_index = index.copy() new_index[masked_locs] = index[unmasked_locs[prev_value]] return new_index def posify_mask_indexer(indexer): """Convert masked values (-1) in an indexer to nearest unmasked values. This routine is useful for dask, where it can be much faster to index adjacent points than arbitrary points from the end of an array. Parameters ---------- indexer : ExplicitIndexer Input indexer. Returns ------- ExplicitIndexer Same type of input, with all values in ndarray keys equal to -1 replaced by an adjacent non-masked element. """ key = tuple( _posify_mask_subindexer(k.ravel()).reshape(k.shape) if isinstance(k, np.ndarray) else k for k in indexer.tuple ) return type(indexer)(key) def is_fancy_indexer(indexer: Any) -> bool: """Return False if indexer is a int, slice, a 1-dimensional list, or a 0 or 1-dimensional ndarray; in all other cases return True """ if isinstance(indexer, (int, slice)): return False if isinstance(indexer, np.ndarray): return indexer.ndim > 1 if isinstance(indexer, list): return bool(indexer) and not isinstance(indexer[0], int) return True class NumpyIndexingAdapter(ExplicitlyIndexedNDArrayMixin): """Wrap a NumPy array to use explicit indexing.""" __slots__ = ("array",) def __init__(self, array): # In NumpyIndexingAdapter we only allow to store bare np.ndarray if not isinstance(array, np.ndarray): raise TypeError( "NumpyIndexingAdapter only wraps np.ndarray. " "Trying to wrap {}".format(type(array)) ) self.array = array def _indexing_array_and_key(self, key): if isinstance(key, OuterIndexer): array = self.array key = _outer_to_numpy_indexer(key, self.array.shape) elif isinstance(key, VectorizedIndexer): array = NumpyVIndexAdapter(self.array) key = key.tuple elif isinstance(key, BasicIndexer): array = self.array # We want 0d slices rather than scalars. This is achieved by # appending an ellipsis (see # https://numpy.org/doc/stable/reference/arrays.indexing.html#detailed-notes). key = key.tuple + (Ellipsis,) else: raise TypeError(f"unexpected key type: {type(key)}") return array, key def transpose(self, order): return self.array.transpose(order) def __getitem__(self, key): array, key = self._indexing_array_and_key(key) return array[key] def __setitem__(self, key, value): array, key = self._indexing_array_and_key(key) try: array[key] = value except ValueError: # More informative exception if read-only view if not array.flags.writeable and not array.flags.owndata: raise ValueError( "Assignment destination is a view. " "Do you want to .copy() array first?" ) else: raise class NdArrayLikeIndexingAdapter(NumpyIndexingAdapter): __slots__ = ("array",) def __init__(self, array): if not hasattr(array, "__array_function__"): raise TypeError( "NdArrayLikeIndexingAdapter must wrap an object that " "implements the __array_function__ protocol" ) self.array = array class ArrayApiIndexingAdapter(ExplicitlyIndexedNDArrayMixin): """Wrap an array API array to use explicit indexing.""" __slots__ = ("array",) def __init__(self, array): if not hasattr(array, "__array_namespace__"): raise TypeError( "ArrayApiIndexingAdapter must wrap an object that " "implements the __array_namespace__ protocol" ) self.array = array def __getitem__(self, key): if isinstance(key, BasicIndexer): return self.array[key.tuple] elif isinstance(key, OuterIndexer): # manual orthogonal indexing (implemented like DaskIndexingAdapter) key = key.tuple value = self.array for axis, subkey in reversed(list(enumerate(key))): value = value[(slice(None),) * axis + (subkey, Ellipsis)] return value else: if isinstance(key, VectorizedIndexer): raise TypeError("Vectorized indexing is not supported") else: raise TypeError(f"Unrecognized indexer: {key}") def __setitem__(self, key, value): if isinstance(key, (BasicIndexer, OuterIndexer)): self.array[key.tuple] = value else: if isinstance(key, VectorizedIndexer): raise TypeError("Vectorized indexing is not supported") else: raise TypeError(f"Unrecognized indexer: {key}") def transpose(self, order): xp = self.array.__array_namespace__() return xp.permute_dims(self.array, order) class DaskIndexingAdapter(ExplicitlyIndexedNDArrayMixin): """Wrap a dask array to support explicit indexing.""" __slots__ = ("array",) def __init__(self, array): """This adapter is created in Variable.__getitem__ in Variable._broadcast_indexes. """ self.array = array def __getitem__(self, key): if not isinstance(key, VectorizedIndexer): # if possible, short-circuit when keys are effectively slice(None) # This preserves dask name and passes lazy array equivalence checks # (see duck_array_ops.lazy_array_equiv) rewritten_indexer = False new_indexer = [] for idim, k in enumerate(key.tuple): if isinstance(k, Iterable) and duck_array_ops.array_equiv( k, np.arange(self.array.shape[idim]) ): new_indexer.append(slice(None)) rewritten_indexer = True else: new_indexer.append(k) if rewritten_indexer: key = type(key)(tuple(new_indexer)) if isinstance(key, BasicIndexer): return self.array[key.tuple] elif isinstance(key, VectorizedIndexer): return self.array.vindex[key.tuple] else: assert isinstance(key, OuterIndexer) key = key.tuple try: return self.array[key] except NotImplementedError: # manual orthogonal indexing. # TODO: port this upstream into dask in a saner way. value = self.array for axis, subkey in reversed(list(enumerate(key))): value = value[(slice(None),) * axis + (subkey,)] return value def __setitem__(self, key, value): if mod_version("dask") >= Version("2021.04.1"): if isinstance(key, BasicIndexer): self.array[key.tuple] = value elif isinstance(key, VectorizedIndexer): self.array.vindex[key.tuple] = value elif isinstance(key, OuterIndexer): num_non_slices = sum( 0 if isinstance(k, slice) else 1 for k in key.tuple ) if num_non_slices > 1: raise NotImplementedError( "xarray can't set arrays with multiple " "array indices to dask yet." ) self.array[key.tuple] = value else: raise TypeError( "This variable's data is stored in a dask array, " "and the installed dask version does not support item " "assignment. To assign to this variable, you must either upgrade dask or" "first load the variable into memory explicitly using the .load() " "method or accessing its .values attribute." ) def transpose(self, order): return self.array.transpose(order) class PandasIndexingAdapter(ExplicitlyIndexedNDArrayMixin): """Wrap a pandas.Index to preserve dtypes and handle explicit indexing.""" __slots__ = ("array", "_dtype") def __init__(self, array: pd.Index, dtype: DTypeLike = None): from xarray.core.indexes import safe_cast_to_index self.array = safe_cast_to_index(array) if dtype is None: self._dtype = get_valid_numpy_dtype(array) else: self._dtype = np.dtype(dtype) @property def dtype(self) -> np.dtype: return self._dtype def __array__(self, dtype: DTypeLike = None) -> np.ndarray: if dtype is None: dtype = self.dtype array = self.array if isinstance(array, pd.PeriodIndex): with suppress(AttributeError): # this might not be public API array = array.astype("object") return np.asarray(array.values, dtype=dtype) @property def shape(self) -> tuple[int, ...]: return (len(self.array),) def _convert_scalar(self, item): if item is pd.NaT: # work around the impossibility of casting NaT with asarray # note: it probably would be better in general to return # pd.Timestamp rather np.than datetime64 but this is easier # (for now) item = np.datetime64("NaT", "ns") elif isinstance(item, timedelta): item = np.timedelta64(getattr(item, "value", item), "ns") elif isinstance(item, pd.Timestamp): # Work around for GH: pydata/xarray#1932 and numpy/numpy#10668 # numpy fails to convert pd.Timestamp to np.datetime64[ns] item = np.asarray(item.to_datetime64()) elif self.dtype != object: item = np.asarray(item, dtype=self.dtype) # as for numpy.ndarray indexing, we always want the result to be # a NumPy array. return to_0d_array(item) def __getitem__( self, indexer ) -> ( PandasIndexingAdapter | NumpyIndexingAdapter | np.ndarray | np.datetime64 | np.timedelta64 ): key = indexer.tuple if isinstance(key, tuple) and len(key) == 1: # unpack key so it can index a pandas.Index object (pandas.Index # objects don't like tuples) (key,) = key if getattr(key, "ndim", 0) > 1: # Return np-array if multidimensional return NumpyIndexingAdapter(np.asarray(self))[indexer] result = self.array[key] if isinstance(result, pd.Index): return type(self)(result, dtype=self.dtype) else: return self._convert_scalar(result) def transpose(self, order) -> pd.Index: return self.array # self.array should be always one-dimensional def __repr__(self) -> str: return f"{type(self).__name__}(array={self.array!r}, dtype={self.dtype!r})" def copy(self, deep: bool = True) -> PandasIndexingAdapter: # Not the same as just writing `self.array.copy(deep=deep)`, as # shallow copies of the underlying numpy.ndarrays become deep ones # upon pickling # >>> len(pickle.dumps((self.array, self.array))) # 4000281 # >>> len(pickle.dumps((self.array, self.array.copy(deep=False)))) # 8000341 array = self.array.copy(deep=True) if deep else self.array return type(self)(array, self._dtype) class PandasMultiIndexingAdapter(PandasIndexingAdapter): """Handles explicit indexing for a pandas.MultiIndex. This allows creating one instance for each multi-index level while preserving indexing efficiency (memoized + might reuse another instance with the same multi-index). """ __slots__ = ("array", "_dtype", "level", "adapter") def __init__( self, array: pd.MultiIndex, dtype: DTypeLike = None, level: str | None = None, ): super().__init__(array, dtype) self.level = level def __array__(self, dtype: DTypeLike = None) -> np.ndarray: if dtype is None: dtype = self.dtype if self.level is not None: return np.asarray( self.array.get_level_values(self.level).values, dtype=dtype ) else: return super().__array__(dtype) def _convert_scalar(self, item): if isinstance(item, tuple) and self.level is not None: idx = tuple(self.array.names).index(self.level) item = item[idx] return super()._convert_scalar(item) def __getitem__(self, indexer): result = super().__getitem__(indexer) if isinstance(result, type(self)): result.level = self.level return result def __repr__(self) -> str: if self.level is None: return super().__repr__() else: props = ( f"(array={self.array!r}, level={self.level!r}, dtype={self.dtype!r})" ) return f"{type(self).__name__}{props}" def _get_array_subset(self) -> np.ndarray: # used to speed-up the repr for big multi-indexes threshold = max(100, OPTIONS["display_values_threshold"] + 2) if self.size > threshold: pos = threshold // 2 indices = np.concatenate([np.arange(0, pos), np.arange(-pos, 0)]) subset = self[OuterIndexer((indices,))] else: subset = self return np.asarray(subset) def _repr_inline_(self, max_width: int) -> str: from xarray.core.formatting import format_array_flat if self.level is None: return "MultiIndex" else: return format_array_flat(self._get_array_subset(), max_width) def _repr_html_(self) -> str: from xarray.core.formatting import short_numpy_repr array_repr = short_numpy_repr(self._get_array_subset()) return f"
{escape(array_repr)}
" def copy(self, deep: bool = True) -> PandasMultiIndexingAdapter: # see PandasIndexingAdapter.copy array = self.array.copy(deep=True) if deep else self.array return type(self)(array, self._dtype, self.level)