Repository: pydata/xarray

Issue:
Inconsistency between Sum of NA's and Mean of NA's: resampling gives 0 or 'NA'
#### Problem description

For datamining with xarray there is always the following issue with the resampling-method.  
If i resample e.g. a daily timeseries over one month and if the data are 'NA' at each day, I get zero as a result. That is annoying considering a timeseries of precipitation. **It is definitely a difference if the monthly precipitation is zero for one month (each day zero precipitation) or the monthly precipitation was not measured due to problems with the device (each day NA)** 

#### Data example

I have a dataset with hourly values for 5 month 'fcut'.
```python
<xarray.Dataset>
Dimensions:       (bnds: 2, time: 3672)
Coordinates:
    rlon          float32 22.06
    rlat          float32 5.06
  * time          (time) datetime64[ns] 2006-05-01 2006-05-01T01:00:00 ...
Dimensions without coordinates: bnds
Data variables:
    rotated_pole  int32 1
    time_bnds     (time, bnds) float64 1.304e+07 1.305e+07 1.305e+07 ...
    TOT_PREC      (time) float64 nan nan nan nan nan nan nan nan nan nan nan ...
Attributes:
```
Doing a resample process gives only zero values for each month.

```
In [10]: fcut.resample(dim='time',freq='M',how='sum')
Out[10]: 
<xarray.Dataset>
Dimensions:       (bnds: 2, time: 5)
Coordinates:
  * time          (time) datetime64[ns] 2006-05-31 2006-06-30 2006-07-31 ...
Dimensions without coordinates: bnds
Data variables:
    rotated_pole  (time) int64 1 1 1 1 1
    time_bnds     (time, bnds) float64 1.07e+10 1.07e+10 1.225e+10 1.225e+10 ...
    TOT_PREC      (time) float64 0.0 0.0 0.0 0.0 0.0
```
**But I expect to have NA for each month, as it is the case for the operator 'mean'**


I know that there is an ongoing discussion about that topic (see for example  https://github.com/pandas-dev/pandas/issues/9422).

For earth science it would be nice to have an option telling xarray what to do in case of a sum over values being all NA. Do you see a chance to have a fast fix for that issue in the model code?

Discussion:
shoyer: Very nice!

> In my implementation, bottleneck is not used when skipna=False.
bottleneck would be advantageous when skipna=True as numpy needs to copy the entire array once, but I think numpy's method is still OK if skipna=False.

I think this is correct -- bottleneck does not speed up non-NaN skipping functions.
fujiisoup: > I think it would make sense to restructure this a little bit to have two well defined layers:

> A module of bottleneck/numpy functions that act on numpy arrays only.
> A module of functions that act on numpy or dask arrays (or these could be moved into duck_array_ops).

Could you explain more detail about this idea?
shoyer: >> A module of bottleneck/numpy functions that act on numpy arrays only.
>>A module of functions that act on numpy or dask arrays (or these could be moved into duck_array_ops).

> Could you explain more detail about this idea?

OK, let me try:

1. On numpy arrays, we use bottleneck eqiuvalents of numpy functions when possible because bottleneck is faster than numpy
2. On dask arrays, we use dask equivalents of numpy functions.
3. We also want to add some extra features on top of what numpy/dask/bottleneck provide, e.g., handling of `min_count`

We could implement this with:
- `nputils.nansum()` is equivalent to `numpy.nansum()` but uses `bottleneck.nansum()` internally instead when possible.
- `duck_array_ops.nansum()` uses `numpy_nansum()` or `dask.array.nansum()`, based upon the type of the inputs.
- `duck_array_ops.sum()` uses `numpy.sum()` or `dask.array.sum()`, based upon the type of the inputs.
- `duck_array_ops.sum_with_mincount()` adds `mincount` and `skipna` support and is used in the `Dataset.sum()` implementation. Its is written using `duck_array_ops.nansum()`, `duck_array_ops.sum()`, `duck_array_ops.where()` and `duck_array_ops.isnull()`.

rpnaut: Hello from me hopefully contributing some needfull things.

At first, I would like to comment that I checked out your code.

I ran the following code example:
``
rpnaut: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!IT SEEMS TO BE A BUG!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

If I do the resample process using the old nomenclature:
` data_aggreg = data["TOT_PREC"].resample(dim="time",how="sum",freq="M") `
it works. 
DO we have a bug in xarray?
rpnaut: Okay. Using the resample-nomenclature I tried to compare the results with your modified code and with the tagged version 0.10.7. I did not saw any nan's in the resampled fields although I expect to have over water only nan's.
fujiisoup: @shoyer , thanks for the details. I think I understood your idea.
This sounds a cleaner solution.
I will update the code again, but it will take some more days (or a week).
fujiisoup: Thanks, @rpnaut .
Actually, I'm changing the code also around `sum`. So it looks my change caused the bug you reported.
I think we do not have a good test coverage around the dataset.reduction.

I will add the test also. Thanks again for the testing :)


fujiisoup: I think this is ready for another review.
rpnaut: Thanks, @fujiisoup .

I have good news and i have bad news. Your min_count argument still seems to work only if using the new syntax for resample, i.e. `data.resample($dim=$freq).sum()`
rpnaut: To wrap it up. Your implementation works for timeseries - data. There is something strange with time-space data, which should be fixed. If this is fixed, it is worth to test in my evaluation environment.
Do you have a feeling, why the new syntax is giving such strange behaviour? Shall we put the bug onto the issue list?

And maybe, it would be interesting to have in the future the min_count argument also available for datasets and not only dataarrays.
fujiisoup: I noticed `min_count` is working also on `resample`d object.
[Your issue](#2356) might be for the API of `resample`.
fujiisoup: Can anyone give further review?

Patch:
diff --git a/.travis.yml b/.travis.yml
index 951b151d..0e51e946 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -95,9 +95,7 @@ install:
   - python xarray/util/print_versions.py
 
 script:
-  # TODO: restore this check once the upstream pandas issue is fixed:
-  # https://github.com/pandas-dev/pandas/issues/21071
-  # - python -OO -c "import xarray"
+  - python -OO -c "import xarray"
   - if [[ "$CONDA_ENV" == "docs" ]]; then
       conda install -c conda-forge sphinx sphinx_rtd_theme sphinx-gallery numpydoc;
       sphinx-build -n -j auto -b html -d _build/doctrees doc _build/html;
diff --git a/asv_bench/benchmarks/unstacking.py b/asv_bench/benchmarks/unstacking.py
new file mode 100644
index 00000000..aa304d4e
--- /dev/null
+++ b/asv_bench/benchmarks/unstacking.py
@@ -0,0 +1,25 @@
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+import xarray as xr
+
+from . import requires_dask
+
+
+class Unstacking(object):
+    def setup(self):
+        data = np.random.RandomState(0).randn(1, 1000, 500)
+        self.ds = xr.DataArray(data).stack(flat_dim=['dim_1', 'dim_2'])
+
+    def time_unstack_fast(self):
+        self.ds.unstack('flat_dim')
+
+    def time_unstack_slow(self):
+        self.ds[:, ::-1].unstack('flat_dim')
+
+
+class UnstackingDask(Unstacking):
+    def setup(self, *args, **kwargs):
+        requires_dask()
+        super(UnstackingDask, self).setup(**kwargs)
+        self.ds = self.ds.chunk({'flat_dim': 50})
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 0c62cca6..4725fe74 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -36,6 +36,12 @@ Documentation
 Enhancements
 ~~~~~~~~~~~~
 
+- min_count option is newly supported in :py:meth:`~xarray.DataArray.sum`,
+  :py:meth:`~xarray.DataArray.prod` and :py:meth:`~xarray.Dataset.sum`, and
+  :py:meth:`~xarray.Dataset.prod`.
+  (:issue:`2230`)
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
+
 - :py:meth:`plot()` now accepts the kwargs ``xscale, yscale, xlim, ylim, xticks, yticks`` just like Pandas. Also ``xincrease=False, yincrease=False`` now use matplotlib's axis inverting methods instead of setting limits.
   By `Deepak Cherian <https://github.com/dcherian>`_. (:issue:`2224`)
 
@@ -52,6 +58,15 @@ Enhancements
   (:issue:`2331`)
   By `Maximilian Roos <https://github.com/max-sixty>`_.
 
+- Applying ``unstack`` to a large DataArray or Dataset is now much faster if the MultiIndex has not been modified after stacking the indices.
+  (:issue:`1560`)
+  By `Maximilian Maahn <https://github.com/maahn>`_.
+
+- You can now control whether or not to offset the coordinates when using
+  the ``roll`` method and the current behavior, coordinates rolled by default,
+  raises a deprecation warning unless explicitly setting the keyword argument.
+  (:issue:`1875`)
+  By `Andrew Huang <https://github.com/ahuang11>`_.
 
 Bug fixes
 ~~~~~~~~~
@@ -61,10 +76,17 @@ Bug fixes
   attribute being set.
   (:issue:`2201`)
   By `Thomas Voigt <https://github.com/tv3141>`_.
+- Fixed a bug in ``zarr`` backend which prevented use with datasets with
+  invalid chunk size encoding after reading from an existing store
+  (:issue:`2278`).
+  By `Joe Hamman <https://github.com/jhamman>`_.
 
 - Tests can be run in parallel with pytest-xdist
   By `Tony Tung <https://github.com/ttung>`_.
 
+- Follow up the renamings in dask; from dask.ghost to dask.overlap 
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
+
 - Now raises a ValueError when there is a conflict between dimension names and
   level names of MultiIndex. (:issue:`2299`)
   By `Keisuke Fujii <https://github.com/fujiisoup>`_.
diff --git a/xarray/backends/api.py b/xarray/backends/api.py
index b2c0df7b..2bf13011 100644
--- a/xarray/backends/api.py
+++ b/xarray/backends/api.py
@@ -806,7 +806,7 @@ def save_mfdataset(datasets, paths, mode='w', format=None, groups=None,
     for obj in datasets:
         if not isinstance(obj, Dataset):
             raise TypeError('save_mfdataset only supports writing Dataset '
-                            'objects, recieved type %s' % type(obj))
+                            'objects, received type %s' % type(obj))
 
     if groups is None:
         groups = [None] * len(datasets)
diff --git a/xarray/backends/zarr.py b/xarray/backends/zarr.py
index a64ca858..47b90c8a 100644
--- a/xarray/backends/zarr.py
+++ b/xarray/backends/zarr.py
@@ -102,9 +102,8 @@ def _determine_zarr_chunks(enc_chunks, var_chunks, ndim):
         enc_chunks_tuple = tuple(enc_chunks)
 
     if len(enc_chunks_tuple) != ndim:
-        raise ValueError("zarr chunks tuple %r must have same length as "
-                         "variable.ndim %g" %
-                         (enc_chunks_tuple, ndim))
+        # throw away encoding chunks, start over
+        return _determine_zarr_chunks(None, var_chunks, ndim)
 
     for x in enc_chunks_tuple:
         if not isinstance(x, int):
diff --git a/xarray/core/common.py b/xarray/core/common.py
index 3f934fcc..55aca5f5 100644
--- a/xarray/core/common.py
+++ b/xarray/core/common.py
@@ -2,6 +2,7 @@ from __future__ import absolute_import, division, print_function
 
 import warnings
 from distutils.version import LooseVersion
+from textwrap import dedent
 
 import numpy as np
 import pandas as pd
@@ -27,20 +28,20 @@ class ImplementsArrayReduce(object):
                                    allow_lazy=True, **kwargs)
         return wrapped_func
 
-    _reduce_extra_args_docstring = \
-        """dim : str or sequence of str, optional
+    _reduce_extra_args_docstring = dedent("""\
+        dim : str or sequence of str, optional
             Dimension(s) over which to apply `{name}`.
         axis : int or sequence of int, optional
             Axis(es) over which to apply `{name}`. Only one of the 'dim'
             and 'axis' arguments can be supplied. If neither are supplied, then
-            `{name}` is calculated over axes."""
+            `{name}` is calculated over axes.""")
 
-    _cum_extra_args_docstring = \
-        """dim : str or sequence of str, optional
+    _cum_extra_args_docstring = dedent("""\
+        dim : str or sequence of str, optional
             Dimension over which to apply `{name}`.
         axis : int or sequence of int, optional
             Axis over which to apply `{name}`. Only one of the 'dim'
-            and 'axis' arguments can be supplied."""
+            and 'axis' arguments can be supplied.""")
 
 
 class ImplementsDatasetReduce(object):
@@ -308,12 +309,12 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         assigned : same type as caller
             A new object with the new coordinates in addition to the existing
             data.
-            
+
         Examples
         --------
-        
+
         Convert longitude coordinates from 0-359 to -180-179:
-        
+
         >>> da = xr.DataArray(np.random.rand(4),
         ...                   coords=[np.array([358, 359, 0, 1])],
         ...                   dims='lon')
@@ -445,11 +446,11 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         grouped : GroupBy
             A `GroupBy` object patterned after `pandas.GroupBy` that can be
             iterated over in the form of `(unique_value, grouped_array)` pairs.
-            
+
         Examples
         --------
         Calculate daily anomalies for daily data:
-        
+
         >>> da = xr.DataArray(np.linspace(0, 1826, num=1827),
         ...                   coords=[pd.date_range('1/1/2000', '31/12/2004',
         ...                           freq='D')],
@@ -465,7 +466,7 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         Coordinates:
           * time       (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03 ...
             dayofyear  (time) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
-        
+
         See Also
         --------
         core.groupby.DataArrayGroupBy
@@ -589,7 +590,7 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
                  closed=None, label=None, base=0, keep_attrs=False, **indexer):
         """Returns a Resample object for performing resampling operations.
 
-        Handles both downsampling and upsampling. If any intervals contain no 
+        Handles both downsampling and upsampling. If any intervals contain no
         values from the original object, they will be given the value ``NaN``.
 
         Parameters
@@ -616,11 +617,11 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         -------
         resampled : same type as caller
             This object resampled.
-            
+
         Examples
         --------
         Downsample monthly time-series data to seasonal data:
-        
+
         >>> da = xr.DataArray(np.linspace(0, 11, num=12),
         ...                   coords=[pd.date_range('15/12/1999',
         ...                           periods=12, freq=pd.DateOffset(months=1))],
@@ -637,13 +638,13 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
           * time     (time) datetime64[ns] 1999-12-01 2000-03-01 2000-06-01 2000-09-01
 
         Upsample monthly time-series data to daily data:
-        
+
         >>> da.resample(time='1D').interpolate('linear')
         <xarray.DataArray (time: 337)>
         array([ 0.      ,  0.032258,  0.064516, ..., 10.935484, 10.967742, 11.      ])
         Coordinates:
           * time     (time) datetime64[ns] 1999-12-15 1999-12-16 1999-12-17 ...
-          
+
         References
         ----------
 
@@ -957,8 +958,8 @@ def contains_cftime_datetimes(var):
                     sample = sample.item()
             return isinstance(sample, cftime_datetime)
         else:
-            return False        
-                    
+            return False
+
 
 def _contains_datetime_like_objects(var):
     """Check if a variable contains datetime like objects (either
diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py
index f215bc47..b1be9944 100644
--- a/xarray/core/dataarray.py
+++ b/xarray/core/dataarray.py
@@ -2015,14 +2015,20 @@ class DataArray(AbstractArray, DataWithCoords):
         variable = self.variable.shift(**shifts)
         return self._replace(variable)
 
-    def roll(self, **shifts):
+    def roll(self, shifts=None, roll_coords=None, **shifts_kwargs):
         """Roll this array by an offset along one or more dimensions.
 
-        Unlike shift, roll rotates all variables, including coordinates. The
-        direction of rotation is consistent with :py:func:`numpy.roll`.
+        Unlike shift, roll may rotate all variables, including coordinates
+        if specified. The direction of rotation is consistent with
+        :py:func:`numpy.roll`.
 
         Parameters
         ----------
+        roll_coords : bool
+            Indicates whether to  roll the coordinates by the offset
+            The current default of roll_coords (None, equivalent to True) is
+            deprecated and will change to False in a future version.
+            Explicitly pass roll_coords to silence the warning.
         **shifts : keyword arguments of the form {dim: offset}
             Integer offset to rotate each of the given dimensions. Positive
             offsets roll to the right; negative offsets roll to the left.
@@ -2046,7 +2052,8 @@ class DataArray(AbstractArray, DataWithCoords):
         Coordinates:
           * x        (x) int64 2 0 1
         """
-        ds = self._to_temp_dataset().roll(**shifts)
+        ds = self._to_temp_dataset().roll(
+            shifts=shifts, roll_coords=roll_coords, **shifts_kwargs)
         return self._from_temp_dataset(ds)
 
     @property
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 4b52178a..37544aca 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -2324,7 +2324,12 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords,
                              'a MultiIndex')
 
         full_idx = pd.MultiIndex.from_product(index.levels, names=index.names)
-        obj = self.reindex(copy=False, **{dim: full_idx})
+
+        # take a shortcut in case the MultiIndex was not modified.
+        if index.equals(full_idx):
+            obj = self
+        else:
+            obj = self.reindex(copy=False, **{dim: full_idx})
 
         new_dim_names = index.names
         new_dim_sizes = [lev.size for lev in index.levels]
@@ -3355,18 +3360,28 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords,
 
         return self._replace_vars_and_dims(variables)
 
-    def roll(self, **shifts):
+    def roll(self, shifts=None, roll_coords=None, **shifts_kwargs):
         """Roll this dataset by an offset along one or more dimensions.
 
-        Unlike shift, roll rotates all variables, including coordinates. The
-        direction of rotation is consistent with :py:func:`numpy.roll`.
+        Unlike shift, roll may rotate all variables, including coordinates
+        if specified. The direction of rotation is consistent with
+        :py:func:`numpy.roll`.
 
         Parameters
         ----------
-        **shifts : keyword arguments of the form {dim: offset}
-            Integer offset to rotate each of the given dimensions. Positive
-            offsets roll to the right; negative offsets roll to the left.
 
+        shifts : dict, optional
+            A dict with keys matching dimensions and values given
+            by integers to rotate each of the given dimensions. Positive
+            offsets roll to the right; negative offsets roll to the left.
+        roll_coords : bool
+            Indicates whether to  roll the coordinates by the offset
+            The current default of roll_coords (None, equivalent to True) is
+            deprecated and will change to False in a future version.
+            Explicitly pass roll_coords to silence the warning.
+        **shifts_kwargs : {dim: offset, ...}, optional
+            The keyword arguments form of ``shifts``.
+            One of shifts or shifts_kwargs must be provided.
         Returns
         -------
         rolled : Dataset
@@ -3389,15 +3404,25 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords,
         Data variables:
             foo      (x) object 'd' 'e' 'a' 'b' 'c'
         """
+        shifts = either_dict_or_kwargs(shifts, shifts_kwargs, 'roll')
         invalid = [k for k in shifts if k not in self.dims]
         if invalid:
             raise ValueError("dimensions %r do not exist" % invalid)
 
+        if roll_coords is None:
+            warnings.warn("roll_coords will be set to False in the future."
+                          " Explicitly set roll_coords to silence warning.",
+                          FutureWarning, stacklevel=2)
+            roll_coords = True
+
+        unrolled_vars = () if roll_coords else self.coords
+
         variables = OrderedDict()
-        for name, var in iteritems(self.variables):
-            var_shifts = dict((k, v) for k, v in shifts.items()
-                              if k in var.dims)
-            variables[name] = var.roll(**var_shifts)
+        for k, v in iteritems(self.variables):
+            if k not in unrolled_vars:
+                variables[k] = v.roll(**shifts)
+            else:
+                variables[k] = v
 
         return self._replace_vars_and_dims(variables)
 
diff --git a/xarray/core/dtypes.py b/xarray/core/dtypes.py
index 7326b936..7ad44472 100644
--- a/xarray/core/dtypes.py
+++ b/xarray/core/dtypes.py
@@ -98,6 +98,9 @@ def maybe_promote(dtype):
     return np.dtype(dtype), fill_value
 
 
+NAT_TYPES = (np.datetime64('NaT'), np.timedelta64('NaT'))
+
+
 def get_fill_value(dtype):
     """Return an appropriate fill value for this dtype.
 
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 3bd10506..17eb310f 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -17,14 +17,6 @@ from . import dask_array_ops, dtypes, npcompat, nputils
 from .nputils import nanfirst, nanlast
 from .pycompat import dask_array_type
 
-try:
-    import bottleneck as bn
-    has_bottleneck = True
-except ImportError:
-    # use numpy methods instead
-    bn = np
-    has_bottleneck = False
-
 try:
     import dask.array as dask_array
     from . import dask_array_compat
@@ -175,7 +167,7 @@ def array_notnull_equiv(arr1, arr2):
 def count(data, axis=None):
     """Count the number of non-NA in this array along the given axis or axes
     """
-    return sum(~isnull(data), axis=axis)
+    return np.sum(~isnull(data), axis=axis)
 
 
 def where(condition, x, y):
@@ -213,159 +205,69 @@ def _ignore_warnings_if(condition):
         yield
 
 
-def _nansum_object(value, axis=None, **kwargs):
-    """ In house nansum for object array """
-    value = fillna(value, 0)
-    return _dask_or_eager_func('sum')(value, axis=axis, **kwargs)
-
-
-def _nan_minmax_object(func, get_fill_value, value, axis=None, **kwargs):
-    """ In house nanmin and nanmax for object array """
-    fill_value = get_fill_value(value.dtype)
-    valid_count = count(value, axis=axis)
-    filled_value = fillna(value, fill_value)
-    data = _dask_or_eager_func(func)(filled_value, axis=axis, **kwargs)
-    if not hasattr(data, 'dtype'):  # scalar case
-        data = dtypes.fill_value(value.dtype) if valid_count == 0 else data
-        return np.array(data, dtype=value.dtype)
-    return where_method(data, valid_count != 0)
-
-
-def _nan_argminmax_object(func, get_fill_value, value, axis=None, **kwargs):
-    """ In house nanargmin, nanargmax for object arrays. Always return integer
-    type """
-    fill_value = get_fill_value(value.dtype)
-    valid_count = count(value, axis=axis)
-    value = fillna(value, fill_value)
-    data = _dask_or_eager_func(func)(value, axis=axis, **kwargs)
-    # dask seems return non-integer type
-    if isinstance(value, dask_array_type):
-        data = data.astype(int)
-
-    if (valid_count == 0).any():
-        raise ValueError('All-NaN slice encountered')
-
-    return np.array(data, dtype=int)
-
-
-def _nanmean_ddof_object(ddof, value, axis=None, **kwargs):
-    """ In house nanmean. ddof argument will be used in _nanvar method """
-    valid_count = count(value, axis=axis)
-    value = fillna(value, 0)
-    # As dtype inference is impossible for object dtype, we assume float
-    # https://github.com/dask/dask/issues/3162
-    dtype = kwargs.pop('dtype', None)
-    if dtype is None and value.dtype.kind == 'O':
-        dtype = value.dtype if value.dtype.kind in ['cf'] else float
-
-    data = _dask_or_eager_func('sum')(value, axis=axis, dtype=dtype, **kwargs)
-    data = data / (valid_count - ddof)
-    return where_method(data, valid_count != 0)
-
-
-def _nanvar_object(value, axis=None, **kwargs):
-    ddof = kwargs.pop('ddof', 0)
-    kwargs_mean = kwargs.copy()
-    kwargs_mean.pop('keepdims', None)
-    value_mean = _nanmean_ddof_object(ddof=0, value=value, axis=axis,
-                                      keepdims=True, **kwargs_mean)
-    squared = (value.astype(value_mean.dtype) - value_mean)**2
-    return _nanmean_ddof_object(ddof, squared, axis=axis, **kwargs)
-
-
-_nan_object_funcs = {
-    'sum': _nansum_object,
-    'min': partial(_nan_minmax_object, 'min', dtypes.get_pos_infinity),
-    'max': partial(_nan_minmax_object, 'max', dtypes.get_neg_infinity),
-    'argmin': partial(_nan_argminmax_object, 'argmin',
-                      dtypes.get_pos_infinity),
-    'argmax': partial(_nan_argminmax_object, 'argmax',
-                      dtypes.get_neg_infinity),
-    'mean': partial(_nanmean_ddof_object, 0),
-    'var': _nanvar_object,
-}
-
-
-def _create_nan_agg_method(name, numeric_only=False, np_compat=False,
-                           no_bottleneck=False, coerce_strings=False):
+def _create_nan_agg_method(name, coerce_strings=False):
+    from . import nanops
+
     def f(values, axis=None, skipna=None, **kwargs):
         if kwargs.pop('out', None) is not None:
             raise TypeError('`out` is not valid for {}'.format(name))
 
-        # If dtype is supplied, we use numpy's method.
-        dtype = kwargs.get('dtype', None)
         values = asarray(values)
 
-        # dask requires dtype argument for object dtype
-        if (values.dtype == 'object' and name in ['sum', ]):
-            kwargs['dtype'] = values.dtype if dtype is None else dtype
-
         if coerce_strings and values.dtype.kind in 'SU':
             values = values.astype(object)
 
+        func = None
         if skipna or (skipna is None and values.dtype.kind in 'cfO'):
-            if values.dtype.kind not in ['u', 'i', 'f', 'c']:
-                func = _nan_object_funcs.get(name, None)
-                using_numpy_nan_func = True
-                if func is None or values.dtype.kind not in 'Ob':
-                    raise NotImplementedError(
-                        'skipna=True not yet implemented for %s with dtype %s'
-                        % (name, values.dtype))
-            else:
-                nanname = 'nan' + name
-                if (isinstance(axis, tuple) or not values.dtype.isnative or
-                        no_bottleneck or (dtype is not None and
-                                          np.dtype(dtype) != values.dtype)):
-                    # bottleneck can't handle multiple axis arguments or
-                    # non-native endianness
-                    if np_compat:
-                        eager_module = npcompat
-                    else:
-                        eager_module = np
-                else:
-                    kwargs.pop('dtype', None)
-                    eager_module = bn
-                func = _dask_or_eager_func(nanname, eager_module)
-                using_numpy_nan_func = (eager_module is np or
-                                        eager_module is npcompat)
+            nanname = 'nan' + name
+            func = getattr(nanops, nanname)
         else:
             func = _dask_or_eager_func(name)
-            using_numpy_nan_func = False
-        with _ignore_warnings_if(using_numpy_nan_func):
-            try:
-                return func(values, axis=axis, **kwargs)
-            except AttributeError:
-                if isinstance(values, dask_array_type):
-                    try:  # dask/dask#3133 dask sometimes needs dtype argument
-                        return func(values, axis=axis, dtype=values.dtype,
-                                    **kwargs)
-                    except AttributeError:
-                        msg = '%s is not yet implemented on dask arrays' % name
-                else:
-                    assert using_numpy_nan_func
-                    msg = ('%s is not available with skipna=False with the '
-                           'installed version of numpy; upgrade to numpy 1.12 '
-                           'or newer to use skipna=True or skipna=None' % name)
-                raise NotImplementedError(msg)
-    f.numeric_only = numeric_only
+
+        try:
+            return func(values, axis=axis, **kwargs)
+        except AttributeError:
+            if isinstance(values, dask_array_type):
+                try:  # dask/dask#3133 dask sometimes needs dtype argument
+                    # if func does not accept dtype, then raises TypeError
+                    return func(values, axis=axis, dtype=values.dtype,
+                                **kwargs)
+                except (AttributeError, TypeError):
+                    msg = '%s is not yet implemented on dask arrays' % name
+            else:
+                msg = ('%s is not available with skipna=False with the '
+                       'installed version of numpy; upgrade to numpy 1.12 '
+                       'or newer to use skipna=True or skipna=None' % name)
+            raise NotImplementedError(msg)
+
     f.__name__ = name
     return f
 
 
+# Attributes `numeric_only`, `available_min_count` is used for docs.
+# See ops.inject_reduce_methods
 argmax = _create_nan_agg_method('argmax', coerce_strings=True)
 argmin = _create_nan_agg_method('argmin', coerce_strings=True)
 max = _create_nan_agg_method('max', coerce_strings=True)
 min = _create_nan_agg_method('min', coerce_strings=True)
-sum = _create_nan_agg_method('sum', numeric_only=True)
-mean = _create_nan_agg_method('mean', numeric_only=True)
-std = _create_nan_agg_method('std', numeric_only=True)
-var = _create_nan_agg_method('var', numeric_only=True)
-median = _create_nan_agg_method('median', numeric_only=True)
-prod = _create_nan_agg_method('prod', numeric_only=True, no_bottleneck=True)
-cumprod_1d = _create_nan_agg_method(
-    'cumprod', numeric_only=True, no_bottleneck=True)
-cumsum_1d = _create_nan_agg_method(
-    'cumsum', numeric_only=True, no_bottleneck=True)
+sum = _create_nan_agg_method('sum')
+sum.numeric_only = True
+sum.available_min_count = True
+mean = _create_nan_agg_method('mean')
+mean.numeric_only = True
+std = _create_nan_agg_method('std')
+std.numeric_only = True
+var = _create_nan_agg_method('var')
+var.numeric_only = True
+median = _create_nan_agg_method('median')
+median.numeric_only = True
+prod = _create_nan_agg_method('prod')
+prod.numeric_only = True
+sum.available_min_count = True
+cumprod_1d = _create_nan_agg_method('cumprod')
+cumprod_1d.numeric_only = True
+cumsum_1d = _create_nan_agg_method('cumsum')
+cumsum_1d.numeric_only = True
 
 
 def _nd_cum_func(cum_func, array, axis, **kwargs):
diff --git a/xarray/core/missing.py b/xarray/core/missing.py
index bec9e2e1..232fa185 100644
--- a/xarray/core/missing.py
+++ b/xarray/core/missing.py
@@ -57,7 +57,7 @@ class NumpyInterpolator(BaseInterpolator):
 
         if self.cons_kwargs:
             raise ValueError(
-                'recieved invalid kwargs: %r' % self.cons_kwargs.keys())
+                'received invalid kwargs: %r' % self.cons_kwargs.keys())
 
         if fill_value is None:
             self._left = np.nan
diff --git a/xarray/core/nanops.py b/xarray/core/nanops.py
new file mode 100644
index 00000000..2309ed96
--- /dev/null
+++ b/xarray/core/nanops.py
@@ -0,0 +1,208 @@
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+
+from . import dtypes
+from .pycompat import dask_array_type
+from . duck_array_ops import (count, isnull, fillna, where_method,
+                              _dask_or_eager_func)
+from . import nputils
+
+try:
+    import dask.array as dask_array
+except ImportError:
+    dask_array = None
+
+
+def _replace_nan(a, val):
+    """
+    replace nan in a by val, and returns the replaced array and the nan
+    position
+    """
+    mask = isnull(a)
+    return where_method(val, mask, a), mask
+
+
+def _maybe_null_out(result, axis, mask, min_count=1):
+    """
+    xarray version of pandas.core.nanops._maybe_null_out
+    """
+    if hasattr(axis, '__len__'):  # if tuple or list
+        raise ValueError('min_count is not available for reduction '
+                         'with more than one dimensions.')
+
+    if axis is not None and getattr(result, 'ndim', False):
+        null_mask = (mask.shape[axis] - mask.sum(axis) - min_count) < 0
+        if null_mask.any():
+            dtype, fill_value = dtypes.maybe_promote(result.dtype)
+            result = result.astype(dtype)
+            result[null_mask] = fill_value
+
+    elif getattr(result, 'dtype', None) not in dtypes.NAT_TYPES:
+        null_mask = mask.size - mask.sum()
+        if null_mask < min_count:
+            result = np.nan
+
+    return result
+
+
+def _nan_argminmax_object(func, fill_value, value, axis=None, **kwargs):
+    """ In house nanargmin, nanargmax for object arrays. Always return integer
+    type
+    """
+    valid_count = count(value, axis=axis)
+    value = fillna(value, fill_value)
+    data = _dask_or_eager_func(func)(value, axis=axis, **kwargs)
+
+    # TODO This will evaluate dask arrays and might be costly.
+    if (valid_count == 0).any():
+        raise ValueError('All-NaN slice encountered')
+
+    return data
+
+
+def _nan_minmax_object(func, fill_value, value, axis=None, **kwargs):
+    """ In house nanmin and nanmax for object array """
+    valid_count = count(value, axis=axis)
+    filled_value = fillna(value, fill_value)
+    data = getattr(np, func)(filled_value, axis=axis, **kwargs)
+    if not hasattr(data, 'dtype'):  # scalar case
+        data = dtypes.fill_value(value.dtype) if valid_count == 0 else data
+        return np.array(data, dtype=value.dtype)
+    return where_method(data, valid_count != 0)
+
+
+def nanmin(a, axis=None, out=None):
+    if a.dtype.kind == 'O':
+        return _nan_minmax_object(
+            'min', dtypes.get_pos_infinity(a.dtype), a, axis)
+
+    module = dask_array if isinstance(a, dask_array_type) else nputils
+    return module.nanmin(a, axis=axis)
+
+
+def nanmax(a, axis=None, out=None):
+    if a.dtype.kind == 'O':
+        return _nan_minmax_object(
+            'max', dtypes.get_neg_infinity(a.dtype), a, axis)
+
+    module = dask_array if isinstance(a, dask_array_type) else nputils
+    return module.nanmax(a, axis=axis)
+
+
+def nanargmin(a, axis=None):
+    fill_value = dtypes.get_pos_infinity(a.dtype)
+    if a.dtype.kind == 'O':
+        return _nan_argminmax_object('argmin', fill_value, a, axis=axis)
+    a, mask = _replace_nan(a, fill_value)
+    if isinstance(a, dask_array_type):
+        res = dask_array.argmin(a, axis=axis)
+    else:
+        res = np.argmin(a, axis=axis)
+
+    if mask is not None:
+        mask = mask.all(axis=axis)
+        if mask.any():
+            raise ValueError("All-NaN slice encountered")
+    return res
+
+
+def nanargmax(a, axis=None):
+    fill_value = dtypes.get_neg_infinity(a.dtype)
+    if a.dtype.kind == 'O':
+        return _nan_argminmax_object('argmax', fill_value, a, axis=axis)
+
+    a, mask = _replace_nan(a, fill_value)
+    if isinstance(a, dask_array_type):
+        res = dask_array.argmax(a, axis=axis)
+    else:
+        res = np.argmax(a, axis=axis)
+
+    if mask is not None:
+        mask = mask.all(axis=axis)
+        if mask.any():
+            raise ValueError("All-NaN slice encountered")
+    return res
+
+
+def nansum(a, axis=None, dtype=None, out=None, min_count=None):
+    a, mask = _replace_nan(a, 0)
+    result = _dask_or_eager_func('sum')(a, axis=axis, dtype=dtype)
+    if min_count is not None:
+        return _maybe_null_out(result, axis, mask, min_count)
+    else:
+        return result
+
+
+def _nanmean_ddof_object(ddof, value, axis=None, **kwargs):
+    """ In house nanmean. ddof argument will be used in _nanvar method """
+    from .duck_array_ops import (count, fillna, _dask_or_eager_func,
+                                 where_method)
+
+    valid_count = count(value, axis=axis)
+    value = fillna(value, 0)
+    # As dtype inference is impossible for object dtype, we assume float
+    # https://github.com/dask/dask/issues/3162
+    dtype = kwargs.pop('dtype', None)
+    if dtype is None and value.dtype.kind == 'O':
+        dtype = value.dtype if value.dtype.kind in ['cf'] else float
+
+    data = _dask_or_eager_func('sum')(value, axis=axis, dtype=dtype, **kwargs)
+    data = data / (valid_count - ddof)
+    return where_method(data, valid_count != 0)
+
+
+def nanmean(a, axis=None, dtype=None, out=None):
+    if a.dtype.kind == 'O':
+        return _nanmean_ddof_object(0, a, axis=axis, dtype=dtype)
+
+    if isinstance(a, dask_array_type):
+        return dask_array.nanmean(a, axis=axis, dtype=dtype)
+
+    return np.nanmean(a, axis=axis, dtype=dtype)
+
+
+def nanmedian(a, axis=None, out=None):
+    return _dask_or_eager_func('nanmedian', eager_module=nputils)(a, axis=axis)
+
+
+def _nanvar_object(value, axis=None, **kwargs):
+    ddof = kwargs.pop('ddof', 0)
+    kwargs_mean = kwargs.copy()
+    kwargs_mean.pop('keepdims', None)
+    value_mean = _nanmean_ddof_object(ddof=0, value=value, axis=axis,
+                                      keepdims=True, **kwargs_mean)
+    squared = (value.astype(value_mean.dtype) - value_mean)**2
+    return _nanmean_ddof_object(ddof, squared, axis=axis, **kwargs)
+
+
+def nanvar(a, axis=None, dtype=None, out=None, ddof=0):
+    if a.dtype.kind == 'O':
+        return _nanvar_object(a, axis=axis, dtype=dtype, ddof=ddof)
+
+    return _dask_or_eager_func('nanvar', eager_module=nputils)(
+        a, axis=axis, dtype=dtype, ddof=ddof)
+
+
+def nanstd(a, axis=None, dtype=None, out=None):
+    return _dask_or_eager_func('nanstd', eager_module=nputils)(
+        a, axis=axis, dtype=dtype)
+
+
+def nanprod(a, axis=None, dtype=None, out=None, min_count=None):
+    a, mask = _replace_nan(a, 1)
+    result = _dask_or_eager_func('nanprod')(a, axis=axis, dtype=dtype, out=out)
+    if min_count is not None:
+        return _maybe_null_out(result, axis, mask, min_count)
+    else:
+        return result
+
+
+def nancumsum(a, axis=None, dtype=None, out=None):
+    return _dask_or_eager_func('nancumsum', eager_module=nputils)(
+        a, axis=axis, dtype=dtype)
+
+
+def nancumprod(a, axis=None, dtype=None, out=None):
+    return _dask_or_eager_func('nancumprod', eager_module=nputils)(
+        a, axis=axis, dtype=dtype)
diff --git a/xarray/core/nputils.py b/xarray/core/nputils.py
index 6df2d34b..a8d596ab 100644
--- a/xarray/core/nputils.py
+++ b/xarray/core/nputils.py
@@ -5,6 +5,14 @@ import warnings
 import numpy as np
 import pandas as pd
 
+try:
+    import bottleneck as bn
+    _USE_BOTTLENECK = True
+except ImportError:
+    # use numpy methods instead
+    bn = np
+    _USE_BOTTLENECK = False
+
 
 def _validate_axis(data, axis):
     ndim = data.ndim
@@ -195,3 +203,36 @@ def _rolling_window(a, window, axis=-1):
     rolling = np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides,
                                               writeable=False)
     return np.swapaxes(rolling, -2, axis)
+
+
+def _create_bottleneck_method(name, npmodule=np):
+    def f(values, axis=None, **kwds):
+        dtype = kwds.get('dtype', None)
+        bn_func = getattr(bn, name, None)
+
+        if (_USE_BOTTLENECK and bn_func is not None and
+                not isinstance(axis, tuple) and
+                values.dtype.kind in 'uifc' and
+                values.dtype.isnative and
+                (dtype is None or np.dtype(dtype) == values.dtype)):
+            # bottleneck does not take care dtype, min_count
+            kwds.pop('dtype', None)
+            result = bn_func(values, axis=axis, **kwds)
+        else:
+            result = getattr(npmodule, name)(values, axis=axis, **kwds)
+
+        return result
+
+    f.__name__ = name
+    return f
+
+
+nanmin = _create_bottleneck_method('nanmin')
+nanmax = _create_bottleneck_method('nanmax')
+nanmean = _create_bottleneck_method('nanmean')
+nanmedian = _create_bottleneck_method('nanmedian')
+nanvar = _create_bottleneck_method('nanvar')
+nanstd = _create_bottleneck_method('nanstd')
+nanprod = _create_bottleneck_method('nanprod')
+nancumsum = _create_bottleneck_method('nancumsum')
+nancumprod = _create_bottleneck_method('nancumprod')
diff --git a/xarray/core/ops.py b/xarray/core/ops.py
index d9e8ceb6..a0dd2212 100644
--- a/xarray/core/ops.py
+++ b/xarray/core/ops.py
@@ -86,7 +86,7 @@ skipna : bool, optional
     If True, skip missing values (as marked by NaN). By default, only
     skips missing values for float dtypes; other dtypes either do not
     have a sentinel missing value (int) or skipna=True has not been
-    implemented (object, datetime64 or timedelta64).
+    implemented (object, datetime64 or timedelta64).{min_count_docs}
 keep_attrs : bool, optional
     If True, the attributes (`attrs`) will be copied from the original
     object to the new one.  If False (default), the new object will be
@@ -102,6 +102,12 @@ reduced : {cls}
     indicated dimension(s) removed.
 """
 
+_MINCOUNT_DOCSTRING = """
+min_count : int, default None
+    The required number of valid values to perform the operation.
+    If fewer than min_count non-NA values are present the result will
+    be NA. New in version 0.10.8: Added with the default being None."""
+
 _ROLLING_REDUCE_DOCSTRING_TEMPLATE = """\
 Reduce this {da_or_ds}'s data windows by applying `{name}` along its dimension.
 
@@ -236,11 +242,15 @@ def inject_reduce_methods(cls):
                [('count', duck_array_ops.count, False)])
     for name, f, include_skipna in methods:
         numeric_only = getattr(f, 'numeric_only', False)
+        available_min_count = getattr(f, 'available_min_count', False)
+        min_count_docs = _MINCOUNT_DOCSTRING if available_min_count else ''
+
         func = cls._reduce_method(f, include_skipna, numeric_only)
         func.__name__ = name
         func.__doc__ = _REDUCE_DOCSTRING_TEMPLATE.format(
             name=name, cls=cls.__name__,
-            extra_args=cls._reduce_extra_args_docstring.format(name=name))
+            extra_args=cls._reduce_extra_args_docstring.format(name=name),
+            min_count_docs=min_count_docs)
         setattr(cls, name, func)
 
 


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index a3d56b90..e6de50b9 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -1404,12 +1404,6 @@ class BaseZarrTest(CFEncodedDataTest):
             with self.roundtrip(ds_chunk4) as actual:
                 self.assertEqual((4,), actual['var1'].encoding['chunks'])
 
-        # specify incompatible encoding
-        ds_chunk4['var1'].encoding.update({'chunks': (5, 5)})
-        with pytest.raises(ValueError) as e_info:
-            with self.roundtrip(ds_chunk4) as actual:
-                pass
-        assert e_info.match('chunks')
 
         # TODO: remove this failure once syncronized overlapping writes are
         # supported by xarray
@@ -1522,6 +1516,21 @@ class BaseZarrTest(CFEncodedDataTest):
             with self.open(store) as actual:
                 assert_identical(original, actual)
 
+    def test_encoding_chunksizes(self):
+        # regression test for GH2278
+        # see also test_encoding_chunksizes_unlimited
+        nx, ny, nt = 4, 4, 5
+        original = xr.Dataset({}, coords={'x': np.arange(nx),
+                                          'y': np.arange(ny),
+                                          't': np.arange(nt)})
+        original['v'] = xr.Variable(('x', 'y', 't'), np.zeros((nx, ny, nt)))
+        original = original.chunk({'t': 1, 'x': 2, 'y': 2})
+
+        with self.roundtrip(original) as ds1:
+            assert_equal(ds1, original)
+            with self.roundtrip(ds1.isel(t=0)) as ds2:
+                assert_equal(ds2, original.isel(t=0))
+
 
 @requires_zarr
 class ZarrDictStoreTest(BaseZarrTest, TestCase):
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index 2950e97c..3619688d 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -3099,9 +3099,24 @@ class TestDataArray(TestCase):
             actual = arr.shift(x=offset)
             assert_identical(expected, actual)
 
-    def test_roll(self):
+    def test_roll_coords(self):
         arr = DataArray([1, 2, 3], coords={'x': range(3)}, dims='x')
-        actual = arr.roll(x=1)
+        actual = arr.roll(x=1, roll_coords=True)
+        expected = DataArray([3, 1, 2], coords=[('x', [2, 0, 1])])
+        assert_identical(expected, actual)
+
+    def test_roll_no_coords(self):
+        arr = DataArray([1, 2, 3], coords={'x': range(3)}, dims='x')
+        actual = arr.roll(x=1, roll_coords=False)
+        expected = DataArray([3, 1, 2], coords=[('x', [0, 1, 2])])
+        assert_identical(expected, actual)
+
+    def test_roll_coords_none(self):
+        arr = DataArray([1, 2, 3], coords={'x': range(3)}, dims='x')
+
+        with pytest.warns(FutureWarning):
+            actual = arr.roll(x=1, roll_coords=None)
+
         expected = DataArray([3, 1, 2], coords=[('x', [2, 0, 1])])
         assert_identical(expected, actual)
 
@@ -3329,7 +3344,9 @@ def test_isin(da):
 def test_rolling_iter(da):
 
     rolling_obj = da.rolling(time=7)
-    rolling_obj_mean = rolling_obj.mean()
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', 'Mean of empty slice')
+        rolling_obj_mean = rolling_obj.mean()
 
     assert len(rolling_obj.window_labels) == len(da['time'])
     assert_identical(rolling_obj.window_labels, da['time'])
@@ -3337,8 +3354,10 @@ def test_rolling_iter(da):
     for i, (label, window_da) in enumerate(rolling_obj):
         assert label == da['time'].isel(time=i)
 
-        actual = rolling_obj_mean.isel(time=i)
-        expected = window_da.mean('time')
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore', 'Mean of empty slice')
+            actual = rolling_obj_mean.isel(time=i)
+            expected = window_da.mean('time')
 
         # TODO add assert_allclose_with_nan, which compares nan position
         # as well as the closeness of the values.
diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py
index 08d71d46..d73632c1 100644
--- a/xarray/tests/test_dataset.py
+++ b/xarray/tests/test_dataset.py
@@ -2113,7 +2113,7 @@ class TestDataset(TestCase):
         with raises_regex(ValueError, 'does not have a MultiIndex'):
             ds.unstack('x')
 
-    def test_stack_unstack(self):
+    def test_stack_unstack_fast(self):
         ds = Dataset({'a': ('x', [0, 1]),
                       'b': (('x', 'y'), [[0, 1], [2, 3]]),
                       'x': [0, 1],
@@ -2124,6 +2124,19 @@ class TestDataset(TestCase):
         actual = ds[['b']].stack(z=['x', 'y']).unstack('z')
         assert actual.identical(ds[['b']])
 
+    def test_stack_unstack_slow(self):
+        ds = Dataset({'a': ('x', [0, 1]),
+                      'b': (('x', 'y'), [[0, 1], [2, 3]]),
+                      'x': [0, 1],
+                      'y': ['a', 'b']})
+        stacked = ds.stack(z=['x', 'y'])
+        actual = stacked.isel(z=slice(None, None, -1)).unstack('z')
+        assert actual.broadcast_equals(ds)
+
+        stacked = ds[['b']].stack(z=['x', 'y'])
+        actual = stacked.isel(z=slice(None, None, -1)).unstack('z')
+        assert actual.identical(ds[['b']])
+
     def test_update(self):
         data = create_test_data(seed=0)
         expected = data.copy()
@@ -2727,6 +2740,20 @@ class TestDataset(TestCase):
             result = actual.reduce(method)
             assert_equal(expected, result)
 
+    def test_resample_min_count(self):
+        times = pd.date_range('2000-01-01', freq='6H', periods=10)
+        ds = Dataset({'foo': (['time', 'x', 'y'], np.random.randn(10, 5, 3)),
+                      'bar': ('time', np.random.randn(10), {'meta': 'data'}),
+                      'time': times})
+        # inject nan
+        ds['foo'] = xr.where(ds['foo'] > 2.0, np.nan, ds['foo'])
+
+        actual = ds.resample(time='1D').sum(min_count=1)
+        expected = xr.concat([
+            ds.isel(time=slice(i * 4, (i + 1) * 4)).sum('time', min_count=1)
+            for i in range(3)], dim=actual['time'])
+        assert_equal(expected, actual)
+
     def test_resample_by_mean_with_keep_attrs(self):
         times = pd.date_range('2000-01-01', freq='6H', periods=10)
         ds = Dataset({'foo': (['time', 'x', 'y'], np.random.randn(10, 5, 3)),
@@ -3365,7 +3392,6 @@ class TestDataset(TestCase):
                                  (('dim2', 'time'), ['dim1', 'dim3']),
                                  ((), ['dim1', 'dim2', 'dim3', 'time'])]:
             actual = data.min(dim=reduct).dims
-            print(reduct, actual, expected)
             self.assertItemsEqual(actual, expected)
 
         assert_equal(data.mean(dim=[]), data)
@@ -3420,7 +3446,6 @@ class TestDataset(TestCase):
                 ('time', ['dim1', 'dim2', 'dim3'])
             ]:
                 actual = getattr(data, cumfunc)(dim=reduct).dims
-                print(reduct, actual, expected)
                 self.assertItemsEqual(actual, expected)
 
     def test_reduce_non_numeric(self):
@@ -3849,18 +3874,42 @@ class TestDataset(TestCase):
         with raises_regex(ValueError, 'dimensions'):
             ds.shift(foo=123)
 
-    def test_roll(self):
+    def test_roll_coords(self):
         coords = {'bar': ('x', list('abc')), 'x': [-4, 3, 2]}
         attrs = {'meta': 'data'}
         ds = Dataset({'foo': ('x', [1, 2, 3])}, coords, attrs)
-        actual = ds.roll(x=1)
+        actual = ds.roll(x=1, roll_coords=True)
 
         ex_coords = {'bar': ('x', list('cab')), 'x': [2, -4, 3]}
         expected = Dataset({'foo': ('x', [3, 1, 2])}, ex_coords, attrs)
         assert_identical(expected, actual)
 
         with raises_regex(ValueError, 'dimensions'):
-            ds.roll(foo=123)
+            ds.roll(foo=123, roll_coords=True)
+
+    def test_roll_no_coords(self):
+        coords = {'bar': ('x', list('abc')), 'x': [-4, 3, 2]}
+        attrs = {'meta': 'data'}
+        ds = Dataset({'foo': ('x', [1, 2, 3])}, coords, attrs)
+        actual = ds.roll(x=1, roll_coords=False)
+
+        expected = Dataset({'foo': ('x', [3, 1, 2])}, coords, attrs)
+        assert_identical(expected, actual)
+
+        with raises_regex(ValueError, 'dimensions'):
+            ds.roll(abc=321, roll_coords=False)
+
+    def test_roll_coords_none(self):
+        coords = {'bar': ('x', list('abc')), 'x': [-4, 3, 2]}
+        attrs = {'meta': 'data'}
+        ds = Dataset({'foo': ('x', [1, 2, 3])}, coords, attrs)
+
+        with pytest.warns(FutureWarning):
+            actual = ds.roll(x=1, roll_coords=None)
+
+        ex_coords = {'bar': ('x', list('cab')), 'x': [2, -4, 3]}
+        expected = Dataset({'foo': ('x', [3, 1, 2])}, ex_coords, attrs)
+        assert_identical(expected, actual)
 
     def test_real_and_imag(self):
         attrs = {'foo': 'bar'}
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index f3f93491..3f32fc49 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -1,12 +1,16 @@
 from __future__ import absolute_import, division, print_function
 
+from distutils.version import LooseVersion
+
 import numpy as np
+import pandas as pd
 import pytest
+from textwrap import dedent
 from numpy import array, nan
 import warnings
 
-from xarray import DataArray, concat
-from xarray.core import duck_array_ops
+from xarray import DataArray, Dataset, concat
+from xarray.core import duck_array_ops, dtypes
 from xarray.core.duck_array_ops import (
     array_notnull_equiv, concatenate, count, first, last, mean, rolling_window,
     stack, where)
@@ -100,7 +104,10 @@ class TestOps(TestCase):
         assert_array_equal(result, np.array([1, 'b'], dtype=object))
 
     def test_all_nan_arrays(self):
-        assert np.isnan(mean([np.nan, np.nan]))
+        with warnings.catch_warnings():
+            warnings.filterwarnings('ignore', 'All-NaN slice')
+            warnings.filterwarnings('ignore', 'Mean of empty slice')
+            assert np.isnan(mean([np.nan, np.nan]))
 
 
 def test_cumsum_1d():
@@ -197,10 +204,15 @@ def construct_dataarray(dim_num, dtype, contains_nan, dask):
         array = rng.choice(['a', 'b', 'c', 'd'], size=shapes)
     else:
         raise ValueError
-    da = DataArray(array, dims=dims, coords={'x': np.arange(16)}, name='da')
 
     if contains_nan:
-        da = da.reindex(x=np.arange(20))
+        inds = rng.choice(range(array.size), int(array.size * 0.2))
+        dtype, fill_value = dtypes.maybe_promote(array.dtype)
+        array = array.astype(dtype)
+        array.flat[inds] = fill_value
+
+    da = DataArray(array, dims=dims, coords={'x': np.arange(16)}, name='da')
+
     if dask and has_dask:
         chunks = {d: 4 for d in dims}
         da = da.chunk(chunks)
@@ -234,10 +246,16 @@ def series_reduce(da, func, dim, **kwargs):
         return concat(da1, dim=d)
 
 
+def assert_dask_array(da, dask):
+    if dask and da.ndim > 0:
+        assert isinstance(da.data, dask_array_type)
+
+
 @pytest.mark.parametrize('dim_num', [1, 2])
 @pytest.mark.parametrize('dtype', [float, int, np.float32, np.bool_])
 @pytest.mark.parametrize('dask', [False, True])
 @pytest.mark.parametrize('func', ['sum', 'min', 'max', 'mean', 'var'])
+# TODO test cumsum, cumprod
 @pytest.mark.parametrize('skipna', [False, True])
 @pytest.mark.parametrize('aggdim', [None, 'x'])
 def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
@@ -251,6 +269,9 @@ def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
     if dask and not has_dask:
         pytest.skip('requires dask')
 
+    if dask and skipna is False and dtype in [np.bool_]:
+        pytest.skip('dask does not compute object-typed array')
+
     rtol = 1e-04 if dtype == np.float32 else 1e-05
 
     da = construct_dataarray(dim_num, dtype, contains_nan=True, dask=dask)
@@ -259,6 +280,7 @@ def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
     # TODO: remove these after resolving
     # https://github.com/dask/dask/issues/3245
     with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', 'Mean of empty slice')
         warnings.filterwarnings('ignore', 'All-NaN slice')
         warnings.filterwarnings('ignore', 'invalid value encountered in')
 
@@ -272,6 +294,7 @@ def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
                     expected = getattr(np, func)(da.values, axis=axis)
 
                 actual = getattr(da, func)(skipna=skipna, dim=aggdim)
+                assert_dask_array(actual, dask)
                 assert np.allclose(actual.values, np.array(expected),
                                    rtol=1.0e-4, equal_nan=True)
             except (TypeError, AttributeError, ZeroDivisionError):
@@ -279,14 +302,21 @@ def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
                 # nanmean for object dtype
                 pass
 
-        # make sure the compatiblility with pandas' results.
         actual = getattr(da, func)(skipna=skipna, dim=aggdim)
+
+        # for dask case, make sure the result is the same for numpy backend
+        expected = getattr(da.compute(), func)(skipna=skipna, dim=aggdim)
+        assert_allclose(actual, expected, rtol=rtol)
+
+        # make sure the compatiblility with pandas' results.
         if func == 'var':
             expected = series_reduce(da, func, skipna=skipna, dim=aggdim,
                                      ddof=0)
             assert_allclose(actual, expected, rtol=rtol)
             # also check ddof!=0 case
             actual = getattr(da, func)(skipna=skipna, dim=aggdim, ddof=5)
+            if dask:
+                assert isinstance(da.data, dask_array_type)
             expected = series_reduce(da, func, skipna=skipna, dim=aggdim,
                                      ddof=5)
             assert_allclose(actual, expected, rtol=rtol)
@@ -297,11 +327,14 @@ def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
         # make sure the dtype argument
         if func not in ['max', 'min']:
             actual = getattr(da, func)(skipna=skipna, dim=aggdim, dtype=float)
+            assert_dask_array(actual, dask)
             assert actual.dtype == float
 
         # without nan
         da = construct_dataarray(dim_num, dtype, contains_nan=False, dask=dask)
         actual = getattr(da, func)(skipna=skipna)
+        if dask:
+            assert isinstance(da.data, dask_array_type)
         expected = getattr(np, 'nan{}'.format(func))(da.values)
         if actual.dtype == object:
             assert actual.values == np.array(expected)
@@ -338,13 +371,6 @@ def test_argmin_max(dim_num, dtype, contains_nan, dask, func, skipna, aggdim):
     with warnings.catch_warnings():
         warnings.filterwarnings('ignore', 'All-NaN slice')
 
-        if aggdim == 'y' and contains_nan and skipna:
-            with pytest.raises(ValueError):
-                actual = da.isel(**{
-                    aggdim: getattr(da, 'arg' + func)(
-                        dim=aggdim, skipna=skipna).compute()})
-            return
-
         actual = da.isel(**{aggdim: getattr(da, 'arg' + func)
                             (dim=aggdim, skipna=skipna).compute()})
         expected = getattr(da, func)(dim=aggdim, skipna=skipna)
@@ -354,6 +380,7 @@ def test_argmin_max(dim_num, dtype, contains_nan, dask, func, skipna, aggdim):
 
 def test_argmin_max_error():
     da = construct_dataarray(2, np.bool_, contains_nan=True, dask=False)
+    da[0] = np.nan
     with pytest.raises(ValueError):
         da.argmin(dim='y')
 
@@ -388,3 +415,122 @@ def test_dask_rolling(axis, window, center):
     with pytest.raises(ValueError):
         rolling_window(dx, axis=axis, window=100, center=center,
                        fill_value=np.nan)
+
+
+@pytest.mark.parametrize('dim_num', [1, 2])
+@pytest.mark.parametrize('dtype', [float, int, np.float32, np.bool_])
+@pytest.mark.parametrize('dask', [False, True])
+@pytest.mark.parametrize('func', ['sum', 'prod'])
+@pytest.mark.parametrize('aggdim', [None, 'x'])
+def test_min_count(dim_num, dtype, dask, func, aggdim):
+    if dask and not has_dask:
+        pytest.skip('requires dask')
+
+    da = construct_dataarray(dim_num, dtype, contains_nan=True, dask=dask)
+    min_count = 3
+
+    actual = getattr(da, func)(dim=aggdim, skipna=True, min_count=min_count)
+
+    if LooseVersion(pd.__version__) >= LooseVersion('0.22.0'):
+        # min_count is only implenented in pandas > 0.22
+        expected = series_reduce(da, func, skipna=True, dim=aggdim,
+                                 min_count=min_count)
+        assert_allclose(actual, expected)
+
+    assert_dask_array(actual, dask)
+
+
+@pytest.mark.parametrize('func', ['sum', 'prod'])
+def test_min_count_dataset(func):
+    da = construct_dataarray(2, dtype=float, contains_nan=True, dask=False)
+    ds = Dataset({'var1': da}, coords={'scalar': 0})
+    actual = getattr(ds, func)(dim='x', skipna=True, min_count=3)['var1']
+    expected = getattr(ds['var1'], func)(dim='x', skipna=True, min_count=3)
+    assert_allclose(actual, expected)
+
+
+@pytest.mark.parametrize('dtype', [float, int, np.float32, np.bool_])
+@pytest.mark.parametrize('dask', [False, True])
+@pytest.mark.parametrize('func', ['sum', 'prod'])
+def test_multiple_dims(dtype, dask, func):
+    if dask and not has_dask:
+        pytest.skip('requires dask')
+    da = construct_dataarray(3, dtype, contains_nan=True, dask=dask)
+
+    actual = getattr(da, func)(('x', 'y'))
+    expected = getattr(getattr(da, func)('x'), func)('y')
+    assert_allclose(actual, expected)
+
+
+def test_docs():
+    # with min_count
+    actual = DataArray.sum.__doc__
+    expected = dedent("""\
+        Reduce this DataArray's data by applying `sum` along some dimension(s).
+
+        Parameters
+        ----------
+        dim : str or sequence of str, optional
+            Dimension(s) over which to apply `sum`.
+        axis : int or sequence of int, optional
+            Axis(es) over which to apply `sum`. Only one of the 'dim'
+            and 'axis' arguments can be supplied. If neither are supplied, then
+            `sum` is calculated over axes.
+        skipna : bool, optional
+            If True, skip missing values (as marked by NaN). By default, only
+            skips missing values for float dtypes; other dtypes either do not
+            have a sentinel missing value (int) or skipna=True has not been
+            implemented (object, datetime64 or timedelta64).
+        min_count : int, default None
+            The required number of valid values to perform the operation.
+            If fewer than min_count non-NA values are present the result will
+            be NA. New in version 0.10.8: Added with the default being None.
+        keep_attrs : bool, optional
+            If True, the attributes (`attrs`) will be copied from the original
+            object to the new one.  If False (default), the new object will be
+            returned without attributes.
+        **kwargs : dict
+            Additional keyword arguments passed on to the appropriate array
+            function for calculating `sum` on this object's data.
+
+        Returns
+        -------
+        reduced : DataArray
+            New DataArray object with `sum` applied to its data and the
+            indicated dimension(s) removed.
+        """)
+    assert actual == expected
+
+    # without min_count
+    actual = DataArray.std.__doc__
+    expected = dedent("""\
+        Reduce this DataArray's data by applying `std` along some dimension(s).
+
+        Parameters
+        ----------
+        dim : str or sequence of str, optional
+            Dimension(s) over which to apply `std`.
+        axis : int or sequence of int, optional
+            Axis(es) over which to apply `std`. Only one of the 'dim'
+            and 'axis' arguments can be supplied. If neither are supplied, then
+            `std` is calculated over axes.
+        skipna : bool, optional
+            If True, skip missing values (as marked by NaN). By default, only
+            skips missing values for float dtypes; other dtypes either do not
+            have a sentinel missing value (int) or skipna=True has not been
+            implemented (object, datetime64 or timedelta64).
+        keep_attrs : bool, optional
+            If True, the attributes (`attrs`) will be copied from the original
+            object to the new one.  If False (default), the new object will be
+            returned without attributes.
+        **kwargs : dict
+            Additional keyword arguments passed on to the appropriate array
+            function for calculating `std` on this object's data.
+
+        Returns
+        -------
+        reduced : DataArray
+            New DataArray object with `std` applied to its data and the
+            indicated dimension(s) removed.
+        """)
+    assert actual == expected
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index cdb578af..3db5e6ad 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -1503,8 +1503,8 @@ class TestVariable(TestCase, VariableSubclassTestCases):
         assert_identical(v.all(dim='x'), Variable([], False))
 
         v = Variable('t', pd.date_range('2000-01-01', periods=3))
-        with pytest.raises(NotImplementedError):
-            v.argmax(skipna=True)
+        assert v.argmax(skipna=True) == 2
+
         assert_identical(
             v.max(), Variable([], pd.Timestamp('2000-01-03')))