diff --git a/xarray/core/_reductions.py b/xarray/core/_reductions.py
index f280f87434c..6b32c60fbaf 100644
--- a/xarray/core/_reductions.py
+++ b/xarray/core/_reductions.py
@@ -73,8 +73,19 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.count()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       int64 5
         """
         return self.reduce(
             duck_array_ops.count,
@@ -133,8 +144,19 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.all()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       bool False
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -193,8 +215,19 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.any()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       bool True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -259,12 +292,27 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.max()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.max(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.max,
@@ -330,12 +378,27 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.min()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.min(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.min,
@@ -405,12 +468,27 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.mean()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 1.8
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.mean(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -487,16 +565,35 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.prod()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 12.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.prod(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.prod(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 12.0
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -574,16 +671,35 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.sum()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 9.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.sum(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.sum(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 9.0
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -658,16 +774,35 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.std()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 0.7483
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.std(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.std(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 0.8367
         """
         return self.reduce(
             duck_array_ops.std,
@@ -742,16 +877,35 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.var()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 0.56
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.var(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.var(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 0.7
         """
         return self.reduce(
             duck_array_ops.var,
@@ -822,12 +976,27 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.median()
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.median(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  ()
+        Data variables:
+            da       float64 nan
         """
         return self.reduce(
             duck_array_ops.median,
@@ -901,8 +1070,15 @@ def count(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.count()
+        <xarray.DataArray ()>
+        array(5)
         """
         return self.reduce(
             duck_array_ops.count,
@@ -959,8 +1135,15 @@ def all(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.all()
+        <xarray.DataArray ()>
+        array(False)
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -1017,8 +1200,15 @@ def any(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.any()
+        <xarray.DataArray ()>
+        array(True)
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -1081,12 +1271,21 @@ def max(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.max()
+        <xarray.DataArray ()>
+        array(3.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.max(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
         """
         return self.reduce(
             duck_array_ops.max,
@@ -1150,12 +1349,21 @@ def min(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.min()
+        <xarray.DataArray ()>
+        array(1.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.min(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
         """
         return self.reduce(
             duck_array_ops.min,
@@ -1223,12 +1431,21 @@ def mean(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.mean()
+        <xarray.DataArray ()>
+        array(1.8)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.mean(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -1303,16 +1520,27 @@ def prod(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.prod()
+        <xarray.DataArray ()>
+        array(12.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.prod(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.prod(skipna=True, min_count=2)
+        <xarray.DataArray ()>
+        array(12.)
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -1388,16 +1616,27 @@ def sum(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.sum()
+        <xarray.DataArray ()>
+        array(9.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.sum(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.sum(skipna=True, min_count=2)
+        <xarray.DataArray ()>
+        array(9.)
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -1470,16 +1709,27 @@ def std(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.std()
+        <xarray.DataArray ()>
+        array(0.74833148)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.std(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.std(skipna=True, ddof=1)
+        <xarray.DataArray ()>
+        array(0.83666003)
         """
         return self.reduce(
             duck_array_ops.std,
@@ -1552,16 +1802,27 @@ def var(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.var()
+        <xarray.DataArray ()>
+        array(0.56)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.var(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.var(skipna=True, ddof=1)
+        <xarray.DataArray ()>
+        array(0.7)
         """
         return self.reduce(
             duck_array_ops.var,
@@ -1630,12 +1891,21 @@ def median(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.median()
+        <xarray.DataArray ()>
+        array(2.)
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.median(skipna=False)
+        <xarray.DataArray ()>
+        array(nan)
         """
         return self.reduce(
             duck_array_ops.median,
@@ -1709,8 +1979,21 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").count()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) int64 1 2 2
         """
         return self.reduce(
             duck_array_ops.count,
@@ -1769,8 +2052,21 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.groupby("labels").all()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) bool False True True
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -1829,8 +2125,21 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.groupby("labels").any()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) bool True True True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -1895,12 +2204,31 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").max()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 2.0 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").max(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 2.0 3.0
         """
         return self.reduce(
             duck_array_ops.max,
@@ -1966,12 +2294,31 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").min()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 2.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").min(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 2.0 1.0
         """
         return self.reduce(
             duck_array_ops.min,
@@ -2041,12 +2388,31 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").mean()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").mean(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 2.0 2.0
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -2123,16 +2489,41 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").prod()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 4.0 3.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").prod(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 4.0 3.0
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.groupby("labels").prod(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 4.0 3.0
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -2210,16 +2601,41 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").sum()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 4.0 4.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").sum(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 4.0 4.0
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.groupby("labels").sum(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 4.0 4.0
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -2294,16 +2710,41 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").std()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 0.0 0.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").std(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 0.0 1.0
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.groupby("labels").std(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 0.0 1.414
         """
         return self.reduce(
             duck_array_ops.std,
@@ -2378,16 +2819,41 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").var()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 0.0 0.0 1.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").var(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 0.0 1.0
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.groupby("labels").var(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 0.0 2.0
         """
         return self.reduce(
             duck_array_ops.var,
@@ -2458,12 +2924,31 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.groupby("labels").median()
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.groupby("labels").median(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (labels: 3)
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
+        Data variables:
+            da       (labels) float64 nan 2.0 2.0
         """
         return self.reduce(
             duck_array_ops.median,
@@ -2538,8 +3023,21 @@ def count(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").count()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) int64 1 3 1
         """
         return self.reduce(
             duck_array_ops.count,
@@ -2598,8 +3096,21 @@ def all(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.resample(time="3M").all()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) bool True True False
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -2658,8 +3169,21 @@ def any(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) bool True True True True True False
 
         >>> ds.resample(time="3M").any()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) bool True True True
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -2724,12 +3248,31 @@ def max(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").max()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 3.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").max(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 3.0 nan
         """
         return self.reduce(
             duck_array_ops.max,
@@ -2795,12 +3338,31 @@ def min(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").min()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 1.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").min(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 1.0 nan
         """
         return self.reduce(
             duck_array_ops.min,
@@ -2870,12 +3432,31 @@ def mean(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").mean()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").mean(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 2.0 nan
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -2952,16 +3533,41 @@ def prod(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").prod()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 6.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").prod(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 6.0 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.resample(time="3M").prod(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 nan 6.0 nan
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -3039,16 +3645,41 @@ def sum(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").sum()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 6.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").sum(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 6.0 nan
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> ds.resample(time="3M").sum(skipna=True, min_count=2)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 nan 6.0 nan
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -3123,16 +3754,41 @@ def std(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").std()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 0.0 0.8165 0.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").std(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 0.0 0.8165 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.resample(time="3M").std(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 nan 1.0 nan
         """
         return self.reduce(
             duck_array_ops.std,
@@ -3207,16 +3863,41 @@ def var(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").var()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 0.0 0.6667 0.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").var(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 0.0 0.6667 nan
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> ds.resample(time="3M").var(skipna=True, ddof=1)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 nan 1.0 nan
         """
         return self.reduce(
             duck_array_ops.var,
@@ -3287,12 +3968,31 @@ def median(
         ... )
         >>> ds = xr.Dataset(dict(da=da))
         >>> ds
+        <xarray.Dataset>
+        Dimensions:  (time: 6)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
+        Data variables:
+            da       (time) float64 1.0 2.0 3.0 1.0 2.0 nan
 
         >>> ds.resample(time="3M").median()
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 2.0 2.0
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> ds.resample(time="3M").median(skipna=False)
+        <xarray.Dataset>
+        Dimensions:  (time: 3)
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
+        Data variables:
+            da       (time) float64 1.0 2.0 nan
         """
         return self.reduce(
             duck_array_ops.median,
@@ -3366,8 +4066,17 @@ def count(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").count()
+        <xarray.DataArray (labels: 3)>
+        array([1, 2, 2])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.count,
@@ -3424,8 +4133,17 @@ def all(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").all()
+        <xarray.DataArray (labels: 3)>
+        array([False,  True,  True])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -3482,8 +4200,17 @@ def any(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").any()
+        <xarray.DataArray (labels: 3)>
+        array([ True,  True,  True])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -3546,12 +4273,25 @@ def max(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").max()
+        <xarray.DataArray (labels: 3)>
+        array([1., 2., 3.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").max(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  2.,  3.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.max,
@@ -3615,12 +4355,25 @@ def min(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").min()
+        <xarray.DataArray (labels: 3)>
+        array([1., 2., 1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").min(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  2.,  1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.min,
@@ -3688,12 +4441,25 @@ def mean(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").mean()
+        <xarray.DataArray (labels: 3)>
+        array([1., 2., 2.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").mean(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  2.,  2.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -3768,16 +4534,33 @@ def prod(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").prod()
+        <xarray.DataArray (labels: 3)>
+        array([1., 4., 3.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").prod(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  4.,  3.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.groupby("labels").prod(skipna=True, min_count=2)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  4.,  3.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -3853,16 +4636,33 @@ def sum(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").sum()
+        <xarray.DataArray (labels: 3)>
+        array([1., 4., 4.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").sum(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  4.,  4.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.groupby("labels").sum(skipna=True, min_count=2)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  4.,  4.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -3935,16 +4735,33 @@ def std(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").std()
+        <xarray.DataArray (labels: 3)>
+        array([0., 0., 1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").std(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  0.,  1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.groupby("labels").std(skipna=True, ddof=1)
+        <xarray.DataArray (labels: 3)>
+        array([       nan, 0.        , 1.41421356])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.std,
@@ -4017,16 +4834,33 @@ def var(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").var()
+        <xarray.DataArray (labels: 3)>
+        array([0., 0., 1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").var(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  0.,  1.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.groupby("labels").var(skipna=True, ddof=1)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  0.,  2.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.var,
@@ -4095,12 +4929,25 @@ def median(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.groupby("labels").median()
+        <xarray.DataArray (labels: 3)>
+        array([1., 2., 2.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.groupby("labels").median(skipna=False)
+        <xarray.DataArray (labels: 3)>
+        array([nan,  2.,  2.])
+        Coordinates:
+          * labels   (labels) object 'a' 'b' 'c'
         """
         return self.reduce(
             duck_array_ops.median,
@@ -4173,8 +5020,17 @@ def count(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").count()
+        <xarray.DataArray (time: 3)>
+        array([1, 3, 1])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.count,
@@ -4231,8 +5087,17 @@ def all(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").all()
+        <xarray.DataArray (time: 3)>
+        array([ True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.array_all,
@@ -4289,8 +5154,17 @@ def any(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ True,  True,  True,  True,  True, False])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").any()
+        <xarray.DataArray (time: 3)>
+        array([ True,  True,  True])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.array_any,
@@ -4353,12 +5227,25 @@ def max(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").max()
+        <xarray.DataArray (time: 3)>
+        array([1., 3., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").max(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  3., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.max,
@@ -4422,12 +5309,25 @@ def min(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").min()
+        <xarray.DataArray (time: 3)>
+        array([1., 1., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").min(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  1., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.min,
@@ -4495,12 +5395,25 @@ def mean(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").mean()
+        <xarray.DataArray (time: 3)>
+        array([1., 2., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").mean(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.mean,
@@ -4575,16 +5488,33 @@ def prod(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").prod()
+        <xarray.DataArray (time: 3)>
+        array([1., 6., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").prod(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  6., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.resample(time="3M").prod(skipna=True, min_count=2)
+        <xarray.DataArray (time: 3)>
+        array([nan,  6., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.prod,
@@ -4660,16 +5590,33 @@ def sum(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").sum()
+        <xarray.DataArray (time: 3)>
+        array([1., 6., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").sum(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  6., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``min_count`` for finer control over when NaNs are ignored.
 
         >>> da.resample(time="3M").sum(skipna=True, min_count=2)
+        <xarray.DataArray (time: 3)>
+        array([nan,  6., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.sum,
@@ -4742,16 +5689,33 @@ def std(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").std()
+        <xarray.DataArray (time: 3)>
+        array([0.        , 0.81649658, 0.        ])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").std(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([0.        , 0.81649658,        nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.resample(time="3M").std(skipna=True, ddof=1)
+        <xarray.DataArray (time: 3)>
+        array([nan,  1., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.std,
@@ -4824,16 +5788,33 @@ def var(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").var()
+        <xarray.DataArray (time: 3)>
+        array([0.        , 0.66666667, 0.        ])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").var(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([0.        , 0.66666667,        nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Specify ``ddof=1`` for an unbiased estimate.
 
         >>> da.resample(time="3M").var(skipna=True, ddof=1)
+        <xarray.DataArray (time: 3)>
+        array([nan,  1., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.var,
@@ -4902,12 +5883,25 @@ def median(
         ...     ),
         ... )
         >>> da
+        <xarray.DataArray (time: 6)>
+        array([ 1.,  2.,  3.,  1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-02-28 ... 2001-06-30
+            labels   (time) <U1 'a' 'b' 'c' 'c' 'b' 'a'
 
         >>> da.resample(time="3M").median()
+        <xarray.DataArray (time: 3)>
+        array([1., 2., 2.])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
 
         Use ``skipna`` to control whether NaNs are ignored.
 
         >>> da.resample(time="3M").median(skipna=False)
+        <xarray.DataArray (time: 3)>
+        array([ 1.,  2., nan])
+        Coordinates:
+          * time     (time) datetime64[ns] 2001-01-31 2001-04-30 2001-07-31
         """
         return self.reduce(
             duck_array_ops.median,
diff --git a/xarray/util/generate_reductions.py b/xarray/util/generate_reductions.py
index ea479ccc9c3..4981efcc7e7 100644
--- a/xarray/util/generate_reductions.py
+++ b/xarray/util/generate_reductions.py
@@ -5,7 +5,7 @@
 Usage:
     python xarray/util/generate_reductions.py > xarray/core/_reductions.py
     pytest --doctest-modules xarray/core/_reductions.py --accept || true
-    pytest --doctest-modules xarray/core/_reductions.py --accept
+    pytest --doctest-modules xarray/core/_reductions.py
 
 This requires [pytest-accept](https://github.com/max-sixty/pytest-accept).
 The second run of pytest is deliberate, since the first will return an error
@@ -25,8 +25,8 @@
 from . import duck_array_ops
 
 if TYPE_CHECKING:
-    from .dataset import Dataset
-    from .dataarray import DataArray'''
+    from .dataarray import DataArray
+    from .dataset import Dataset'''
 
 
 CLASS_PREAMBLE = """