Merge branch 'main' into groupby-aggs-using-numpy-groupies

* main: (23 commits)
  Vectorize groupby binary ops (pydata#6160)
  Speed-up multi-index html repr + add display_values_threshold option (pydata#6400)
  [pre-commit.ci] pre-commit autoupdate (pydata#6422)
  Fix concat scalar coord dtype (pydata#6418)
  use the `DaskIndexingAdapter` for `duck dask` arrays (pydata#6414)
  Weighted quantile (pydata#6059)
  upgrade `sphinx` (pydata#6415)
  Add kwarg-only breaking change to whats-new (pydata#6409)
  [pre-commit.ci] pre-commit autoupdate (pydata#6396)
  fix DataArray groupby returning a Dataset (pydata#6394)
  reindex: fix missing variable metadata (pydata#6389)
  [skip-ci] Add benchmarks for groupby math (pydata#6390)
  Fix concat with scalar coordinate (pydata#6385)
  isel: convert IndexVariable to Variable if index is dropped (pydata#6388)
  fix dataset groupby combine dataarray func (pydata#6386)
  fix concat with variable or dataarray as dim (pydata#6387)
  pydata#6367 Fix for time units checking could produce "unhashable type" error (pydata#6368)
  Explicit indexes (pydata#5692)
  Remove test_rasterio_vrt_network (pydata#6371)
  Allow write_empty_chunks to be set in Zarr encoding (pydata#6348)
  ...

Author: dcherian

.pre-commit-config.yaml

@@ -19,7 +19,7 @@ repos:
     hooks:
       - id: isort
   - repo: https://github.com/asottile/pyupgrade
-    rev: v2.31.0
+    rev: v2.31.1
     hooks:
       - id: pyupgrade
         args:
@@ -45,7 +45,7 @@ repos:
   #     - id: velin
   #       args: ["--write", "--compact"]
   - repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v0.931
+    rev: v0.942
     hooks:
       - id: mypy
         # Copied from setup.cfg

readthedocs.yml renamed to .readthedocs.yaml

@@ -1,10 +1,14 @@
 version: 2
+
 build:
   os: ubuntu-20.04
   tools:
     python: mambaforge-4.10
-sphinx:
-  fail_on_warning: true
+
 conda:
   environment: ci/requirements/doc.yml
+
+sphinx:
+  fail_on_warning: true
+
 formats: []

asv_bench/benchmarks/groupby.py

@@ -17,6 +17,8 @@ def setup(self, *args, **kwargs):
             }
         )
         self.ds2d = self.ds1d.expand_dims(z=10)
+        self.ds1d_mean = self.ds1d.groupby("b").mean()
+        self.ds2d_mean = self.ds2d.groupby("b").mean()
 
     @parameterized(["ndim"], [(1, 2)])
     def time_init(self, ndim):
@@ -32,15 +34,30 @@ def time_agg_large_num_groups(self, method, ndim):
         ds = getattr(self, f"ds{ndim}d")
         getattr(ds.groupby("b"), method)()
 
+    def time_groupby_binary_op_1d(self):
+        self.ds1d - self.ds1d_mean
+
+    def time_groupby_binary_op_2d(self):
+        self.ds2d - self.ds2d_mean
+
+    def peakmem_groupby_binary_op_1d(self):
+        self.ds1d - self.ds1d_mean
+
+    def peakmem_groupby_binary_op_2d(self):
+        self.ds2d - self.ds2d_mean
+
 
 class GroupByDask(GroupBy):
     def setup(self, *args, **kwargs):
         requires_dask()
         super().setup(**kwargs)
+
         self.ds1d = self.ds1d.sel(dim_0=slice(None, None, 2))
         self.ds1d["c"] = self.ds1d["c"].chunk({"dim_0": 50})
         self.ds2d = self.ds2d.sel(dim_0=slice(None, None, 2))
         self.ds2d["c"] = self.ds2d["c"].chunk({"dim_0": 50, "z": 5})
+        self.ds1d_mean = self.ds1d.groupby("b").mean()
+        self.ds2d_mean = self.ds2d.groupby("b").mean()
 
 
 class GroupByPandasDataFrame(GroupBy):
@@ -52,6 +69,13 @@ def setup(self, *args, **kwargs):
 
         super().setup(**kwargs)
         self.ds1d = self.ds1d.to_dataframe()
+        self.ds1d_mean = self.ds1d.groupby("b").mean()
+
+    def time_groupby_binary_op_2d(self):
+        raise NotImplementedError
+
+    def peakmem_groupby_binary_op_2d(self):
+        raise NotImplementedError
 
 
 class GroupByDaskDataFrame(GroupBy):
@@ -64,6 +88,13 @@ def setup(self, *args, **kwargs):
         requires_dask()
         super().setup(**kwargs)
         self.ds1d = self.ds1d.chunk({"dim_0": 50}).to_dataframe()
+        self.ds1d_mean = self.ds1d.groupby("b").mean()
+
+    def time_groupby_binary_op_2d(self):
+        raise NotImplementedError
+
+    def peakmem_groupby_binary_op_2d(self):
+        raise NotImplementedError
 
 
 class Resample:
@@ -75,6 +106,8 @@ def setup(self, *args, **kwargs):
             coords={"time": pd.date_range("2001-01-01", freq="H", periods=365 * 24)},
         )
         self.ds2d = self.ds1d.expand_dims(z=10)
+        self.ds1d_mean = self.ds1d.resample(time="48H").mean()
+        self.ds2d_mean = self.ds2d.resample(time="48H").mean()
 
     @parameterized(["ndim"], [(1, 2)])
     def time_init(self, ndim):
@@ -90,6 +123,18 @@ def time_agg_large_num_groups(self, method, ndim):
         ds = getattr(self, f"ds{ndim}d")
         getattr(ds.resample(time="48H"), method)()
 
+    def time_groupby_binary_op_1d(self):
+        self.ds1d - self.ds1d_mean
+
+    def time_groupby_binary_op_2d(self):
+        self.ds2d - self.ds2d_mean
+
+    def peakmem_groupby_binary_op_1d(self):
+        self.ds1d - self.ds1d_mean
+
+    def peakmem_groupby_binary_op_2d(self):
+        self.ds2d - self.ds2d_mean
+
 
 class ResampleDask(Resample):
     def setup(self, *args, **kwargs):

ci/requirements/doc.yml

@@ -4,16 +4,16 @@ channels:
   - conda-forge
   - nodefaults
 dependencies:
-  - python=3.8
+  - python=3.9
   - bottleneck
   - cartopy
   - cfgrib>=0.9
   - dask-core>=2.30
   - h5netcdf>=0.7.4
   - ipykernel
   - ipython
-  - ipython_genutils  # remove once `nbconvert` fixed its dependencies
   - iris>=2.3
+  - jinja2<3.1  # remove once nbconvert fixed the use of removed functions
   - jupyter_client
   - matplotlib-base
   - nbsphinx
@@ -34,7 +34,7 @@ dependencies:
   - sphinx-book-theme >= 0.0.38
   - sphinx-copybutton
   - sphinx-panels
-  - sphinx<4
+  - sphinx!=4.4.0
   - zarr>=2.4
   - pip:
       - sphinxext-rediraffe

doc/api.rst

@@ -944,6 +944,7 @@ Dataset
 
    DatasetWeighted
    DatasetWeighted.mean
+   DatasetWeighted.quantile
    DatasetWeighted.sum
    DatasetWeighted.std
    DatasetWeighted.var
@@ -958,6 +959,7 @@ DataArray
 
    DataArrayWeighted
    DataArrayWeighted.mean
+   DataArrayWeighted.quantile
    DataArrayWeighted.sum
    DataArrayWeighted.std
    DataArrayWeighted.var

doc/developers-meeting.rst

@@ -7,7 +7,7 @@ The meeting occurs on `Zoom <https://us02web.zoom.us/j/88251613296?pwd=azZsSkU1U
 
 Notes for the meeting are kept `here <https://hackmd.io/@U4W-olO3TX-hc-cvbjNe4A/xarray-dev-meeting/edit>`__.
 
-There is a `GitHub issue <https://github.com/pydata/xarray/issues/4001>`__ for changes to the meeting.
+There is a :issue:`GitHub issue <4001>` for changes to the meeting.
 
 You can subscribe to this calendar to be notified of changes:
 

doc/internals/extending-xarray.rst

@@ -18,12 +18,10 @@ easy to inadvertently use internal APIs when subclassing, which means that your
 code may break when xarray upgrades. Furthermore, many builtin methods will
 only return native xarray objects.
 
-The standard advice is to use `composition over inheritance`__, but
+The standard advice is to use :issue:`composition over inheritance <706>`, but
 reimplementing an API as large as xarray's on your own objects can be an onerous
 task, even if most methods are only forwarding to xarray implementations.
 
-__ https://github.com/pydata/xarray/issues/706
-
 If you simply want the ability to call a function with the syntax of a
 method call, then the builtin :py:meth:`~xarray.DataArray.pipe` method (copied
 from pandas) may suffice.

doc/roadmap.rst

@@ -114,8 +114,7 @@ xarray's data model, e.g., as attributes on the ``Dataset`` and
 coordinates in xarray operations, but will no longer would need to have
 a one-to-one correspondence with coordinate variables. Instead, an index
 should be able to refer to multiple (possibly multidimensional)
-coordinates that define it. See `GH
-1603 <https://github.com/pydata/xarray/issues/1603>`__ for full details
+coordinates that define it. See :issue:`1603` for full details.
 
 Specific tasks:
 
@@ -182,11 +181,9 @@ backends means that users can not easily build backend interface for
 xarray in third-party libraries.
 
 The idea of refactoring the backends API and exposing it to users was
-originally proposed in `GH
-1970 <https://github.com/pydata/xarray/issues/1970>`__. The idea would
-be to develop a well tested and generic backend base class and
-associated utilities for external use. Specific tasks for this
-development would include:
+originally proposed in :issue:`1970`. The idea would be to develop a
+well tested and generic backend base class and associated utilities
+for external use. Specific tasks for this development would include:
 
 -  Exposing an abstract backend for writing new storage systems.
 -  Exposing utilities for features like automatic closing of files,
@@ -225,7 +222,7 @@ examples include:
 
 A new tree-like data structure which is essentially a structured hierarchical
 collection of Datasets could represent these cases, and would instead map to
-multiple netCDF groups (see `GH4118 <https://github.com/pydata/xarray/issues/4118>`__.).
+multiple netCDF groups (see :issue:`4118`).
 
 Currently there are several libraries which have wrapped xarray in order to build
 domain-specific data structures (e.g. `xarray-multiscale <https://github.com/JaneliaSciComp/xarray-multiscale>`__.),

doc/tutorials-and-videos.rst

@@ -18,6 +18,20 @@ Videos
 .. panels::
   :card: text-center
 
+  ---
+  Xdev Python Tutorial Seminar Series 2022 Thinking with Xarray : High-level computation patterns | Deepak Cherian
+  ^^^
+  .. raw:: html
+
+    <iframe width="100%" src="https://www.youtube.com/embed/TSw3GF_d2y8" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
+
+  ---
+  Xdev Python Tutorial Seminar Series 2021 seminar introducing xarray (2 of 2) | Anderson Banihirwe
+  ^^^
+  .. raw:: html
+
+    <iframe width="100%" src="https://www.youtube.com/embed/2H_4drBwORY" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
+
   ---
   Xdev Python Tutorial Seminar Series 2021 seminar introducing xarray (1 of 2) | Anderson Banihirwe
   ^^^

doc/user-guide/computation.rst

@@ -265,7 +265,7 @@ Weighted array reductions
 
 :py:class:`DataArray` and :py:class:`Dataset` objects include :py:meth:`DataArray.weighted`
 and :py:meth:`Dataset.weighted` array reduction methods. They currently
-support weighted ``sum``, ``mean``, ``std`` and ``var``.
+support weighted ``sum``, ``mean``, ``std``, ``var`` and ``quantile``.
 
 .. ipython:: python
 
@@ -293,6 +293,12 @@ Calculate the weighted mean:
 
     weighted_prec.mean(dim="month")
 
+Calculate the weighted quantile:
+
+.. ipython:: python
+
+    weighted_prec.quantile(q=0.5, dim="month")
+
 The weighted sum corresponds to:
 
 .. ipython:: python

doc/user-guide/plotting.rst

@@ -251,7 +251,7 @@ Finally, if a dataset does not have any coordinates it enumerates all data point
 .. ipython:: python
     :okwarning:
 
-    air1d_multi = air1d_multi.drop("date")
+    air1d_multi = air1d_multi.drop(["date", "time", "decimal_day"])
     air1d_multi.plot()
 
 The same applies to 2D plots below.

doc/whats-new.rst

@@ -22,10 +22,30 @@ v2022.03.1 (unreleased)
 New Features
 ~~~~~~~~~~~~
 
+- Add a weighted ``quantile`` method to :py:class:`~core.weighted.DatasetWeighted` and
+  :py:class:`~core.weighted.DataArrayWeighted` (:pull:`6059`). By
+  `Christian Jauvin <https://github.com/cjauvin>`_ and `David Huard <https://github.com/huard>`_.
+- Add a ``create_index=True`` parameter to :py:meth:`Dataset.stack` and
+  :py:meth:`DataArray.stack` so that the creation of multi-indexes is optional
+  (:pull:`5692`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
+- Multi-index levels are now accessible through their own, regular coordinates
+  instead of virtual coordinates (:pull:`5692`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
+- Add a ``display_values_threshold`` option to control the total number of array
+  elements which trigger summarization rather than full repr in (numpy) array
+  detailed views of the html repr (:pull:`6400`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
 
+- The Dataset and DataArray ``rename*`` methods do not implicitly add or drop
+  indexes. (:pull:`5692`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
+- Many arguments like ``keep_attrs``, ``axis``, and ``skipna`` are now keyword
+  only for all reduction operations like ``.mean``.
+  By `Deepak Cherian <https://github.com/dcherian>`_, `Jimmy Westling <https://github.com/illviljan>`_.
 
 Deprecations
 ~~~~~~~~~~~~
@@ -36,17 +56,36 @@ Bug fixes
 
 - Set ``skipna=None`` for all ``quantile`` methods (e.g. :py:meth:`Dataset.quantile`) and
   ensure it skips missing values for float dtypes (consistent with other methods). This should
-  not change the behavior (:pull:`6303`). By `Mathias Hauser <https://github.com/mathause>`_.
+  not change the behavior (:pull:`6303`).
+  By `Mathias Hauser <https://github.com/mathause>`_.
+- Many bugs fixed by the explicit indexes refactor, mainly related to multi-index (virtual)
+  coordinates. See the corresponding pull-request on GitHub for more details. (:pull:`5692`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
+- Fixed "unhashable type" error trying to read NetCDF file with variable having its 'units'
+  attribute not ``str`` (e.g. ``numpy.ndarray``) (:issue:`6368`).
+  By `Oleh Khoma <https://github.com/okhoma>`_.
+- Fixed the poor html repr performance on large multi-indexes (:pull:`6400`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
+- Allow fancy indexing of duck dask arrays along multiple dimensions. (:pull:`6414`)
+  By `Justus Magin <https://github.com/keewis>`_.
 
 Documentation
 ~~~~~~~~~~~~~
 
+Performance
+~~~~~~~~~~~
+
+- GroupBy binary operations are now vectorized.
+  Previously this involved looping over all groups. (:issue:`5804`,:pull:`6160`)
+  By `Deepak Cherian <https://github.com/dcherian>`_.
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
 
+- Many internal changes due to the explicit indexes refactor. See the
+  corresponding pull-request on GitHub for more details. (:pull:`5692`).
+  By `Benoît Bovy <https://github.com/benbovy>`_.
 
-.. _whats-new.2022.02.0:
 .. _whats-new.2022.03.0:
 
 v2022.03.0 (2 March 2022)
@@ -2129,7 +2168,7 @@ Documentation
 - Created a "How do I..." section (:ref:`howdoi`) for solutions to common questions. (:pull:`3357`).
   By `Deepak Cherian <https://github.com/dcherian>`_.
 - Add examples for :py:meth:`Dataset.swap_dims` and :py:meth:`DataArray.swap_dims`
-  (pull:`3331`, pull:`3331`). By `Justus Magin <https://github.com/keewis>`_.
+  (:pull:`3331`, :pull:`3331`). By `Justus Magin <https://github.com/keewis>`_.
 - Add examples for :py:meth:`align`, :py:meth:`merge`, :py:meth:`combine_by_coords`,
   :py:meth:`full_like`, :py:meth:`zeros_like`, :py:meth:`ones_like`, :py:meth:`Dataset.pipe`,
   :py:meth:`Dataset.assign`, :py:meth:`Dataset.reindex`, :py:meth:`Dataset.fillna` (:pull:`3328`).
@@ -2713,7 +2752,7 @@ Removes inadvertently introduced setup dependency on pytest-runner
   will be Python 3 only, but older versions of xarray will always be available
   for Python 2.7 users. For the more details, see:
 
-  - `Xarray Github issue discussing dropping Python 2 <https://github.com/pydata/xarray/issues/1829>`__
+  - :issue:`Xarray Github issue discussing dropping Python 2 <1829>`
   - `Python 3 Statement <http://www.python3statement.org/>`__
   - `Tips on porting to Python 3 <https://docs.python.org/3/howto/pyporting.html>`__
 

xarray/backends/zarr.py

@@ -212,7 +212,13 @@ def extract_zarr_variable_encoding(
     """
     encoding = variable.encoding.copy()
 
-    valid_encodings = {"chunks", "compressor", "filters", "cache_metadata"}
+    valid_encodings = {
+        "chunks",
+        "compressor",
+        "filters",
+        "cache_metadata",
+        "write_empty_chunks",
+    }
 
     if raise_on_invalid:
         invalid = [k for k in encoding if k not in valid_encodings]

xarray/coding/times.py

@@ -695,7 +695,8 @@ def encode(self, variable, name=None):
     def decode(self, variable, name=None):
         dims, data, attrs, encoding = unpack_for_decoding(variable)
 
-        if "units" in attrs and attrs["units"] in TIME_UNITS:
+        units = attrs.get("units")
+        if isinstance(units, str) and units in TIME_UNITS:
             units = pop_to(attrs, encoding, "units")
             transform = partial(decode_cf_timedelta, units=units)
             dtype = np.dtype("timedelta64[ns]")