Repository: pydata/xarray

Issue:
Zarr keys include variable name
When using open_zarr on a dataset with many variables the keynames include the variable name, like 

    ('zarr-temperature-1234', 1, 3, 2)

In the distributed scheduler these keynames get shortened to prefixes like `zarr-temperature` and used both for scheduling heuristics (all keys with the same prefix are expected to take similar-ish amounts of time) and for diagnostics, such as in the progress bar plot below:

![image](https://user-images.githubusercontent.com/306380/35920009-9f1fbeda-0be4-11e8-8f3f-061387bdb149.png)

We may want to avoid including the variable name into the keyname here in order to avoid breaking these out into several groups.  Instead you might consider putting the variable name within the key as another member of the tuple like the following:

    ('zarr-1234', 'temperature', 1, 3 ,2)

Discussion:

Patch:
diff --git a/.stickler.yml b/.stickler.yml
new file mode 100644
index 00000000..1634c1f5
--- /dev/null
+++ b/.stickler.yml
@@ -0,0 +1,7 @@
+linters:
+  flake8:
+    max-line-length: 79
+    fixer: true
+  py3k:
+fixers:
+  enable: true
diff --git a/.travis.yml b/.travis.yml
index 068ea3cc..ee8ffcc4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -45,6 +45,8 @@ matrix:
     env: CONDA_ENV=py36-rasterio1.0alpha
   - python: 3.6
     env: CONDA_ENV=py36-zarr-dev
+  - python: 3.5
+    env: CONDA_ENV=docs
   allow_failures:
   - python: 3.6
     env:
@@ -86,16 +88,25 @@ before_install:
   - conda info -a
 
 install:
-  - conda env create --file ci/requirements-$CONDA_ENV.yml
+  - if [[ "$CONDA_ENV" == "docs" ]]; then
+      conda env create -n test_env --file doc/environment.yml;
+    else
+      conda env create -n test_env --file ci/requirements-$CONDA_ENV.yml;
+    fi
   - source activate test_env
   - conda list
-  - python setup.py install
+  - pip install --no-deps -e .
   - python xarray/util/print_versions.py
 
 script:
   - flake8 -j auto xarray
   - python -OO -c "import xarray"
-  - py.test xarray --cov=xarray --cov-config ci/.coveragerc --cov-report term-missing --verbose $EXTRA_FLAGS
+  - if [[ "$CONDA_ENV" == "docs" ]]; then
+      conda install -c conda-forge sphinx_rtd_theme;
+      sphinx-build -n -b html -d _build/doctrees doc _build/html;
+    else
+      py.test xarray --cov=xarray --cov-config ci/.coveragerc --cov-report term-missing --verbose $EXTRA_FLAGS;
+    fi
 
 after_success:
   - coveralls
diff --git a/asv_bench/asv.conf.json b/asv_bench/asv.conf.json
index a2878a7b..b5953436 100644
--- a/asv_bench/asv.conf.json
+++ b/asv_bench/asv.conf.json
@@ -63,7 +63,7 @@
         "netcdf4": [""],
         "scipy": [""],
         "bottleneck": ["", null],
-        "dask": ["", null],
+        "dask": [""],
     },
 
 
diff --git a/asv_bench/benchmarks/reindexing.py b/asv_bench/benchmarks/reindexing.py
new file mode 100644
index 00000000..0f28eaa4
--- /dev/null
+++ b/asv_bench/benchmarks/reindexing.py
@@ -0,0 +1,45 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+
+import numpy as np
+import xarray as xr
+
+from . import requires_dask
+
+
+class Reindex(object):
+    def setup(self):
+        data = np.random.RandomState(0).randn(1000, 100, 100)
+        self.ds = xr.Dataset({'temperature': (('time', 'x', 'y'), data)},
+                             coords={'time': np.arange(1000),
+                                     'x': np.arange(100),
+                                     'y': np.arange(100)})
+
+    def time_1d_coarse(self):
+        self.ds.reindex(time=np.arange(0, 1000, 5)).load()
+
+    def time_1d_fine_all_found(self):
+        self.ds.reindex(time=np.arange(0, 1000, 0.5), method='nearest').load()
+
+    def time_1d_fine_some_missing(self):
+        self.ds.reindex(time=np.arange(0, 1000, 0.5), method='nearest',
+                        tolerance=0.1).load()
+
+    def time_2d_coarse(self):
+        self.ds.reindex(x=np.arange(0, 100, 2), y=np.arange(0, 100, 2)).load()
+
+    def time_2d_fine_all_found(self):
+        self.ds.reindex(x=np.arange(0, 100, 0.5), y=np.arange(0, 100, 0.5),
+                        method='nearest').load()
+
+    def time_2d_fine_some_missing(self):
+        self.ds.reindex(x=np.arange(0, 100, 0.5), y=np.arange(0, 100, 0.5),
+                        method='nearest', tolerance=0.1).load()
+
+
+class ReindexDask(Reindex):
+    def setup(self):
+        requires_dask()
+        super(ReindexDask, self).setup()
+        self.ds = self.ds.chunk({'time': 100})
diff --git a/doc/environment.yml b/doc/environment.yml
index b14fba35..2758612c 100644
--- a/doc/environment.yml
+++ b/doc/environment.yml
@@ -18,3 +18,4 @@ dependencies:
   - sphinx-gallery
   - zarr
   - iris
+  - flake8
diff --git a/doc/installing.rst b/doc/installing.rst
index b6ec2bd8..b9a1fff5 100644
--- a/doc/installing.rst
+++ b/doc/installing.rst
@@ -78,6 +78,7 @@ Testing
 -------
 
 To run the test suite after installing xarray, first install (via pypi or conda)
+
 - `py.test <https://pytest.org>`__: Simple unit testing library
 - `mock <https://pypi.python.org/pypi/mock>`__: additional testing library required for python version 2
 
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 875aba51..fe8332c3 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -43,7 +43,16 @@ Documentation
 
 Enhancements
 ~~~~~~~~~~~~
-- reduce methods such as :py:func:`DataArray.sum()` now accepts ``dtype``
+- Reduce methods such as :py:func:`DataArray.sum()` now handles object-type array.
+
+  .. ipython:: python
+
+    da = xr.DataArray(np.array([True, False, np.nan], dtype=object), dims='x')
+    da.sum()
+
+  (:issue:`1866`)
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
+- Reduce methods such as :py:func:`DataArray.sum()` now accepts ``dtype``
   arguments. (:issue:`1838`)
   By `Keisuke Fujii <https://github.com/fujiisoup>`_.
 - Added nodatavals attribute to DataArray when using :py:func:`~xarray.open_rasterio`. (:issue:`1736`).
@@ -81,6 +90,9 @@ Enhancements
   By `Zac Hatfield-Dodds <https://github.com/Zac-HD>`_.
 - `.dt` accessor can now ceil, floor and round timestamps to specified frequency.
   By `Deepak Cherian <https://github.com/dcherian>`_.
+- Speed of reindexing/alignment with dask array is orders of magnitude faster
+  when inserting missing values  (:issue:`1847`).
+  By `Stephan Hoyer <https://github.com/shoyer>`_.
 
 .. _Zarr: http://zarr.readthedocs.io/
 
@@ -140,6 +152,13 @@ Bug fixes
   ``parse_coordinates`` kwarg has beed added to :py:func:`~open_rasterio`
   (set to ``True`` per default).
   By `Fabien Maussion <https://github.com/fmaussion>`_.
+- The colors of discrete colormaps are now the same regardless if `seaborn`
+  is installed or not (:issue:`1896`).
+  By `Fabien Maussion <https://github.com/fmaussion>`_.
+- Fixed dtype promotion rules in :py:func:`where` and :py:func:`concat` to
+  match pandas (:issue:`1847`). A combination of strings/numbers or
+  unicode/bytes now promote to object dtype, instead of strings or unicode.
+  By `Stephan Hoyer <https://github.com/shoyer>`_.
 
 .. _whats-new.0.10.0:
 
@@ -476,6 +495,9 @@ Bug fixes
   ``apionly`` module was deprecated.
   (:issue:`1633`). By `Joe Hamman <https://github.com/jhamman>`_.
 
+- Fix COMPAT: MultiIndex checking is fragile
+  (:issue:`1833`). By `Florian Pinault <https://github.com/floriankrb>`_.
+
 - Fix ``rasterio`` backend for Rasterio versions 1.0alpha10 and newer.
   (:issue:`1641`). By `Chris Holden <https://github.com/ceholden>`_.
 
diff --git a/xarray/backends/zarr.py b/xarray/backends/zarr.py
index 02753f6c..2737d9fb 100644
--- a/xarray/backends/zarr.py
+++ b/xarray/backends/zarr.py
@@ -482,7 +482,7 @@ def open_zarr(store, group=None, synchronizer=None, auto_chunk=True,
             if (var.ndim > 0) and (chunks is not None):
                 # does this cause any data to be read?
                 token2 = tokenize(name, var._data)
-                name2 = 'zarr-%s-%s' % (name, token2)
+                name2 = 'zarr-%s' % token2
                 return var.chunk(chunks, name=name2, lock=None)
             else:
                 return var
diff --git a/xarray/core/alignment.py b/xarray/core/alignment.py
index 87624532..99dde45b 100644
--- a/xarray/core/alignment.py
+++ b/xarray/core/alignment.py
@@ -8,13 +8,11 @@ import warnings
 
 import numpy as np
 
-from . import duck_array_ops
-from . import dtypes
 from . import utils
 from .indexing import get_indexer_nd
 from .pycompat import iteritems, OrderedDict, suppress
 from .utils import is_full_slice, is_dict_like
-from .variable import Variable, IndexVariable
+from .variable import IndexVariable
 
 
 def _get_joiner(join):
@@ -306,59 +304,51 @@ def reindex_variables(variables, sizes, indexes, indexers, method=None,
     from .dataarray import DataArray
 
     # build up indexers for assignment along each dimension
-    to_indexers = {}
-    from_indexers = {}
+    int_indexers = {}
+    targets = {}
+    masked_dims = set()
+    unchanged_dims = set()
+
     # size of reindexed dimensions
     new_sizes = {}
 
     for name, index in iteritems(indexes):
         if name in indexers:
-            target = utils.safe_cast_to_index(indexers[name])
             if not index.is_unique:
                 raise ValueError(
                     'cannot reindex or align along dimension %r because the '
                     'index has duplicate values' % name)
-            indexer = get_indexer_nd(index, target, method, tolerance)
 
+            target = utils.safe_cast_to_index(indexers[name])
             new_sizes[name] = len(target)
-            # Note pandas uses negative values from get_indexer_nd to signify
-            # values that are missing in the index
-            # The non-negative values thus indicate the non-missing values
-            to_indexers[name] = indexer >= 0
-            if to_indexers[name].all():
-                # If an indexer includes no negative values, then the
-                # assignment can be to a full-slice (which is much faster,
-                # and means we won't need to fill in any missing values)
-                to_indexers[name] = slice(None)
-
-            from_indexers[name] = indexer[to_indexers[name]]
-            if np.array_equal(from_indexers[name], np.arange(len(index))):
-                # If the indexer is equal to the original index, use a full
-                # slice object to speed up selection and so we can avoid
-                # unnecessary copies
-                from_indexers[name] = slice(None)
+
+            int_indexer = get_indexer_nd(index, target, method, tolerance)
+
+            # We uses negative values from get_indexer_nd to signify
+            # values that are missing in the index.
+            if (int_indexer < 0).any():
+                masked_dims.add(name)
+            elif np.array_equal(int_indexer, np.arange(len(index))):
+                unchanged_dims.add(name)
+
+            int_indexers[name] = int_indexer
+            targets[name] = target
 
     for dim in sizes:
         if dim not in indexes and dim in indexers:
             existing_size = sizes[dim]
-            new_size = utils.safe_cast_to_index(indexers[dim]).size
+            new_size = indexers[dim].size
             if existing_size != new_size:
                 raise ValueError(
                     'cannot reindex or align along dimension %r without an '
                     'index because its size %r is different from the size of '
                     'the new index %r' % (dim, existing_size, new_size))
 
-    def any_not_full_slices(indexers):
-        return any(not is_full_slice(idx) for idx in indexers)
-
-    def var_indexers(var, indexers):
-        return tuple(indexers.get(d, slice(None)) for d in var.dims)
-
     # create variables for the new dataset
     reindexed = OrderedDict()
 
     for dim, indexer in indexers.items():
-        if isinstance(indexer, DataArray) and indexer.dims != (dim, ):
+        if isinstance(indexer, DataArray) and indexer.dims != (dim,):
             warnings.warn(
                 "Indexer has dimensions {0:s} that are different "
                 "from that to be indexed along {1:s}. "
@@ -375,47 +365,24 @@ def reindex_variables(variables, sizes, indexes, indexers, method=None,
 
     for name, var in iteritems(variables):
         if name not in indexers:
-            assign_to = var_indexers(var, to_indexers)
-            assign_from = var_indexers(var, from_indexers)
-
-            if any_not_full_slices(assign_to):
-                # there are missing values to in-fill
-                data = var[assign_from].data
-                dtype, fill_value = dtypes.maybe_promote(var.dtype)
-
-                if isinstance(data, np.ndarray):
-                    shape = tuple(new_sizes.get(dim, size)
-                                  for dim, size in zip(var.dims, var.shape))
-                    new_data = np.empty(shape, dtype=dtype)
-                    new_data[...] = fill_value
-                    # create a new Variable so we can use orthogonal indexing
-                    # use fastpath=True to avoid dtype inference
-                    new_var = Variable(var.dims, new_data, var.attrs,
-                                       fastpath=True)
-                    new_var[assign_to] = data
-
-                else:  # dask array
-                    data = data.astype(dtype, copy=False)
-                    for axis, indexer in enumerate(assign_to):
-                        if not is_full_slice(indexer):
-                            indices = np.cumsum(indexer)[~indexer]
-                            data = duck_array_ops.insert(
-                                data, indices, fill_value, axis=axis)
-                    new_var = Variable(var.dims, data, var.attrs,
-                                       fastpath=True)
-
-            elif any_not_full_slices(assign_from):
-                # type coercion is not necessary as there are no missing
-                # values
-                new_var = var[assign_from]
-
-            else:
-                # no reindexing is necessary
+            key = tuple(slice(None)
+                        if d in unchanged_dims
+                        else int_indexers.get(d, slice(None))
+                        for d in var.dims)
+            needs_masking = any(d in masked_dims for d in var.dims)
+
+            if needs_masking:
+                new_var = var._getitem_with_mask(key)
+            elif all(is_full_slice(k) for k in key):
+                # no reindexing necessary
                 # here we need to manually deal with copying data, since
                 # we neither created a new ndarray nor used fancy indexing
                 new_var = var.copy(deep=copy)
+            else:
+                new_var = var[key]
 
             reindexed[name] = new_var
+
     return reindexed
 
 
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 62ad2b9b..5f31a3d9 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -1253,7 +1253,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, BaseDataObject,
             from dask.base import tokenize
         except ImportError:
             import dask  # raise the usual error if dask is entirely missing  # flake8: noqa
-            raise ImportError('xarray requires dask version 0.6 or newer')
+            raise ImportError('xarray requires dask version 0.9 or newer')
 
         if isinstance(chunks, Number):
             chunks = dict.fromkeys(self.dims, chunks)
@@ -2752,7 +2752,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, BaseDataObject,
         idx = dataframe.index
         obj = cls()
 
-        if hasattr(idx, 'levels'):
+        if isinstance(idx, pd.MultiIndex):
             # it's a multi-index
             # expand the DataFrame to include the product of all levels
             full_idx = pd.MultiIndex.from_product(idx.levels, names=idx.names)
diff --git a/xarray/core/dtypes.py b/xarray/core/dtypes.py
index ccbe48ed..8dac3961 100644
--- a/xarray/core/dtypes.py
+++ b/xarray/core/dtypes.py
@@ -1,4 +1,5 @@
 import numpy as np
+import functools
 
 from . import utils
 
@@ -7,6 +8,63 @@ from . import utils
 NA = utils.ReprObject('<NA>')
 
 
+@functools.total_ordering
+class AlwaysGreaterThan(object):
+    def __gt__(self, other):
+        return True
+
+    def __eq__(self, other):
+        return isinstance(other, type(self))
+
+
+@functools.total_ordering
+class AlwaysLessThan(object):
+    def __lt__(self, other):
+        return True
+
+    def __eq__(self, other):
+        return isinstance(other, type(self))
+
+
+# Equivalence to np.inf (-np.inf) for object-type
+INF = AlwaysGreaterThan()
+NINF = AlwaysLessThan()
+
+
+# Pairs of types that, if both found, should be promoted to object dtype
+# instead of following NumPy's own type-promotion rules. These type promotion
+# rules match pandas instead. For reference, see the NumPy type hierarchy:
+# https://docs.scipy.org/doc/numpy-1.13.0/reference/arrays.scalars.html
+PROMOTE_TO_OBJECT = [
+    {np.number, np.character},  # numpy promotes to character
+    {np.bool_, np.character},  # numpy promotes to character
+    {np.bytes_, np.unicode_},  # numpy promotes to unicode
+]
+
+
+@functools.total_ordering
+class AlwaysGreaterThan(object):
+    def __gt__(self, other):
+        return True
+
+    def __eq__(self, other):
+        return isinstance(other, type(self))
+
+
+@functools.total_ordering
+class AlwaysLessThan(object):
+    def __lt__(self, other):
+        return True
+
+    def __eq__(self, other):
+        return isinstance(other, type(self))
+
+
+# Equivalence to np.inf (-np.inf) for object-type
+INF = AlwaysGreaterThan()
+NINF = AlwaysLessThan()
+
+
 def maybe_promote(dtype):
     """Simpler equivalent of pandas.core.common._maybe_promote
 
@@ -55,8 +113,74 @@ def get_fill_value(dtype):
     return fill_value
 
 
+def get_pos_infinity(dtype):
+    """Return an appropriate positive infinity for this dtype.
+
+    Parameters
+    ----------
+    dtype : np.dtype
+
+    Returns
+    -------
+    fill_value : positive infinity value corresponding to this dtype.
+    """
+    if issubclass(dtype.type, (np.floating, np.integer)):
+        return np.inf
+
+    if issubclass(dtype.type, np.complexfloating):
+        return np.inf + 1j * np.inf
+
+    return INF
+
+
+def get_neg_infinity(dtype):
+    """Return an appropriate positive infinity for this dtype.
+
+    Parameters
+    ----------
+    dtype : np.dtype
+
+    Returns
+    -------
+    fill_value : positive infinity value corresponding to this dtype.
+    """
+    if issubclass(dtype.type, (np.floating, np.integer)):
+        return -np.inf
+
+    if issubclass(dtype.type, np.complexfloating):
+        return -np.inf - 1j * np.inf
+
+    return NINF
+
+
 def is_datetime_like(dtype):
     """Check if a dtype is a subclass of the numpy datetime types
     """
     return (np.issubdtype(dtype, np.datetime64) or
             np.issubdtype(dtype, np.timedelta64))
+
+
+def result_type(*arrays_and_dtypes):
+    """Like np.result_type, but with type promotion rules matching pandas.
+
+    Examples of changed behavior:
+    number + string -> object (not string)
+    bytes + unicode -> object (not unicode)
+
+    Parameters
+    ----------
+    *arrays_and_dtypes : list of arrays and dtypes
+        The dtype is extracted from both numpy and dask arrays.
+
+    Returns
+    -------
+    numpy.dtype for the result.
+    """
+    types = {np.result_type(t).type for t in arrays_and_dtypes}
+
+    for left, right in PROMOTE_TO_OBJECT:
+        if (any(issubclass(t, left) for t in types) and
+                any(issubclass(t, right) for t in types)):
+            return np.dtype(object)
+
+    return np.result_type(*arrays_and_dtypes)
diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py
index 2058ce86..6f554880 100644
--- a/xarray/core/duck_array_ops.py
+++ b/xarray/core/duck_array_ops.py
@@ -82,13 +82,13 @@ def isnull(data):
 
 
 transpose = _dask_or_eager_func('transpose')
-where = _dask_or_eager_func('where', n_array_args=3)
+_where = _dask_or_eager_func('where', n_array_args=3)
 insert = _dask_or_eager_func('insert')
 take = _dask_or_eager_func('take')
 broadcast_to = _dask_or_eager_func('broadcast_to')
 
-concatenate = _dask_or_eager_func('concatenate', list_of_args=True)
-stack = _dask_or_eager_func('stack', list_of_args=True)
+_concatenate = _dask_or_eager_func('concatenate', list_of_args=True)
+_stack = _dask_or_eager_func('stack', list_of_args=True)
 
 array_all = _dask_or_eager_func('all')
 array_any = _dask_or_eager_func('any')
@@ -100,6 +100,17 @@ def asarray(data):
     return data if isinstance(data, dask_array_type) else np.asarray(data)
 
 
+def as_shared_dtype(scalars_or_arrays):
+    """Cast a arrays to a shared dtype using xarray's type promotion rules."""
+    arrays = [asarray(x) for x in scalars_or_arrays]
+    # Pass arrays directly instead of dtypes to result_type so scalars
+    # get handled properly.
+    # Note that result_type() safely gets the dtype from dask arrays without
+    # evaluating them.
+    out_type = dtypes.result_type(*arrays)
+    return [x.astype(out_type, copy=False) for x in arrays]
+
+
 def as_like_arrays(*data):
     if all(isinstance(d, dask_array_type) for d in data):
         return data
@@ -151,6 +162,11 @@ def count(data, axis=None):
     return sum(~isnull(data), axis=axis)
 
 
+def where(condition, x, y):
+    """Three argument where() with better dtype promotion rules."""
+    return _where(condition, *as_shared_dtype([x, y]))
+
+
 def where_method(data, cond, other=dtypes.NA):
     if other is dtypes.NA:
         other = dtypes.get_fill_value(data.dtype)
@@ -161,6 +177,16 @@ def fillna(data, other):
     return where(isnull(data), other, data)
 
 
+def concatenate(arrays, axis=0):
+    """concatenate() with better dtype promotion rules."""
+    return _concatenate(as_shared_dtype(arrays), axis=axis)
+
+
+def stack(arrays, axis=0):
+    """stack() with better dtype promotion rules."""
+    return _stack(as_shared_dtype(arrays), axis=axis)
+
+
 @contextlib.contextmanager
 def _ignore_warnings_if(condition):
     if condition:
@@ -171,6 +197,79 @@ 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,
                            keep_dims=False):
@@ -185,27 +284,31 @@ def _create_nan_agg_method(name, numeric_only=False, np_compat=False,
         if coerce_strings and values.dtype.kind in 'SU':
             values = values.astype(object)
 
-        if skipna or (skipna is None and values.dtype.kind in 'cf'):
+        if skipna or (skipna is None and values.dtype.kind in 'cfO'):
             if values.dtype.kind not in ['u', 'i', 'f', 'c']:
-                raise NotImplementedError(
-                    'skipna=True not yet implemented for %s with dtype %s'
-                    % (name, values.dtype))
-            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
+                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:
-                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
+                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)
         else:
             func = _dask_or_eager_func(name)
             using_numpy_nan_func = False
@@ -214,7 +317,11 @@ def _create_nan_agg_method(name, numeric_only=False, np_compat=False,
                 return func(values, axis=axis, **kwargs)
             except AttributeError:
                 if isinstance(values, dask_array_type):
-                    msg = '%s is not yet implemented on dask arrays' % name
+                    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 '
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index d4863014..14fc3480 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -1273,8 +1273,6 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
 
         arrays = [v.data for v in variables]
 
-        # TODO: use our own type promotion rules to ensure that
-        # [str, float] -> object, not str like numpy
         if dim in first_var.dims:
             axis = first_var.get_axis_num(dim)
             dims = first_var.dims
diff --git a/xarray/plot/utils.py b/xarray/plot/utils.py
index c194b9dd..0e565f24 100644
--- a/xarray/plot/utils.py
+++ b/xarray/plot/utils.py
@@ -115,32 +115,23 @@ def _color_palette(cmap, n_colors):
     colors_i = np.linspace(0, 1., n_colors)
     if isinstance(cmap, (list, tuple)):
         # we have a list of colors
-        try:
-            sns = import_seaborn()
-        except ImportError:
-            # if that fails, use matplotlib
-            # in this case, is there any difference between mpl and seaborn?
-            cmap = ListedColormap(cmap, N=n_colors)
-            pal = cmap(colors_i)
-        else:
-            # first try to turn it into a palette with seaborn
-            pal = sns.color_palette(cmap, n_colors=n_colors)
+        cmap = ListedColormap(cmap, N=n_colors)
+        pal = cmap(colors_i)
     elif isinstance(cmap, basestring):
         # we have some sort of named palette
         try:
-            # first try to turn it into a palette with seaborn
-            from seaborn.apionly import color_palette
-            pal = color_palette(cmap, n_colors=n_colors)
-        except (ImportError, ValueError):
-            # ValueError is raised when seaborn doesn't like a colormap
-            # (e.g. jet). If that fails, use matplotlib
+            # is this a matplotlib cmap?
+            cmap = plt.get_cmap(cmap)
+            pal = cmap(colors_i)
+        except ValueError:
+            # ValueError happens when mpl doesn't like a colormap, try seaborn
             try:
-                # is this a matplotlib cmap?
-                cmap = plt.get_cmap(cmap)
-            except ValueError:
+                from seaborn.apionly import color_palette
+                pal = color_palette(cmap, n_colors=n_colors)
+            except (ValueError, ImportError):
                 # or maybe we just got a single color as a string
                 cmap = ListedColormap([cmap], N=n_colors)
-            pal = cmap(colors_i)
+                pal = cmap(colors_i)
     else:
         # cmap better be a LinearSegmentedColormap (e.g. viridis)
         pal = cmap(colors_i)


Test patch:
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index 39b3109c..0def5b68 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -1717,6 +1717,12 @@ class TestDataArray(TestCase):
         actual = arr.where(arr.x < 2, drop=True)
         assert_identical(actual, expected)
 
+    def test_where_string(self):
+        array = DataArray(['a', 'b'])
+        expected = DataArray(np.array(['a', np.nan], dtype=object))
+        actual = array.where([True, False])
+        assert_identical(actual, expected)
+
     def test_cumops(self):
         coords = {'x': [-1, -2], 'y': ['ab', 'cd', 'ef'],
                   'lat': (['x', 'y'], [[1, 2, 3], [-1, -2, -3]]),
diff --git a/xarray/tests/test_dtypes.py b/xarray/tests/test_dtypes.py
new file mode 100644
index 00000000..1b236e01
--- /dev/null
+++ b/xarray/tests/test_dtypes.py
@@ -0,0 +1,54 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+
+import numpy as np
+import pytest
+
+from xarray.core import dtypes
+
+
+@pytest.mark.parametrize("args, expected", [
+    ([np.bool], np.bool),
+    ([np.bool, np.string_], np.object_),
+    ([np.float32, np.float64], np.float64),
+    ([np.float32, np.string_], np.object_),
+    ([np.unicode_, np.int64], np.object_),
+    ([np.unicode_, np.unicode_], np.unicode_),
+    ([np.bytes_, np.unicode_], np.object_),
+])
+def test_result_type(args, expected):
+    actual = dtypes.result_type(*args)
+    assert actual == expected
+
+
+def test_result_type_scalar():
+    actual = dtypes.result_type(np.arange(3, dtype=np.float32), np.nan)
+    assert actual == np.float32
+
+
+def test_result_type_dask_array():
+    # verify it works without evaluating dask arrays
+    da = pytest.importorskip('dask.array')
+    dask = pytest.importorskip('dask')
+
+    def error():
+        raise RuntimeError
+
+    array = da.from_delayed(dask.delayed(error)(), (), np.float64)
+    with pytest.raises(RuntimeError):
+        array.compute()
+
+    actual = dtypes.result_type(array)
+    assert actual == np.float64
+
+    # note that this differs from the behavior for scalar numpy arrays, which
+    # would get promoted to float32
+    actual = dtypes.result_type(array, np.array([0.5, 1.0], dtype=np.float32))
+    assert actual == np.float64
+
+
+@pytest.mark.parametrize('obj', [1.0, np.inf, 'ab', 1.0 + 1.0j, True])
+def test_inf(obj):
+    assert dtypes.INF > obj
+    assert dtypes.NINF < obj
diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py
index 9fb1b1aa..d68a7a38 100644
--- a/xarray/tests/test_duck_array_ops.py
+++ b/xarray/tests/test_duck_array_ops.py
@@ -1,15 +1,19 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from pytest import mark
+import pytest
 import numpy as np
 from numpy import array, nan
+from distutils.version import LooseVersion
 from . import assert_array_equal
 from xarray.core.duck_array_ops import (
-    first, last, count, mean, array_notnull_equiv,
+    first, last, count, mean, array_notnull_equiv, where, stack, concatenate
 )
+from xarray import DataArray
+from xarray.testing import assert_allclose
+from xarray import concat
 
-from . import TestCase, raises_regex
+from . import TestCase, raises_regex, has_dask
 
 
 class TestOps(TestCase):
@@ -76,12 +80,28 @@ class TestOps(TestCase):
         expected = array([[1, 2, 3], [3, 2, 1]])
         assert_array_equal(expected, count(self.x, axis=-1))
 
+    def test_where_type_promotion(self):
+        result = where([True, False], [1, 2], ['a', 'b'])
+        assert_array_equal(result, np.array([1, 'b'], dtype=object))
+
+        result = where([True, False], np.array([1, 2], np.float32), np.nan)
+        assert result.dtype == np.float32
+        assert_array_equal(result, np.array([1, np.nan], dtype=np.float32))
+
+    def test_stack_type_promotion(self):
+        result = stack([1, 'b'])
+        assert_array_equal(result, np.array([1, 'b'], dtype=object))
+
+    def test_concatenate_type_promotion(self):
+        result = concatenate([[1], ['b']])
+        assert_array_equal(result, np.array([1, 'b'], dtype=object))
+
     def test_all_nan_arrays(self):
         assert np.isnan(mean([np.nan, np.nan]))
 
 
 class TestArrayNotNullEquiv():
-    @mark.parametrize("arr1, arr2", [
+    @pytest.mark.parametrize("arr1, arr2", [
         (np.array([1, 2, 3]), np.array([1, 2, 3])),
         (np.array([1, 2, np.nan]), np.array([1, np.nan, 3])),
         (np.array([np.nan, 2, np.nan]), np.array([1, np.nan, np.nan])),
@@ -99,7 +119,7 @@ class TestArrayNotNullEquiv():
         b = np.array([[1, 2], [np.nan, 4]])
         assert not array_notnull_equiv(a, b)
 
-    @mark.parametrize("val1, val2, val3, null", [
+    @pytest.mark.parametrize("val1, val2, val3, null", [
         (1, 2, 3, None),
         (1., 2., 3., np.nan),
         (1., 2., 3., None),
@@ -109,3 +129,175 @@ class TestArrayNotNullEquiv():
         arr1 = np.array([val1, null, val3, null])
         arr2 = np.array([val1, val2, null, null])
         assert array_notnull_equiv(arr1, arr2)
+
+
+def construct_dataarray(dim_num, dtype, contains_nan, dask):
+    # dimnum <= 3
+    rng = np.random.RandomState(0)
+    shapes = [16, 8, 4][:dim_num]
+    dims = ('x', 'y', 'z')[:dim_num]
+
+    if np.issubdtype(dtype, np.floating):
+        array = rng.randn(*shapes).astype(dtype)
+    elif np.issubdtype(dtype, np.integer):
+        array = rng.randint(0, 10, size=shapes).astype(dtype)
+    elif np.issubdtype(dtype, np.bool_):
+        array = rng.randint(0, 1, size=shapes).astype(dtype)
+    elif dtype == str:
+        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))
+    if dask and has_dask:
+        chunks = {d: 4 for d in dims}
+        da = da.chunk(chunks)
+
+    return da
+
+
+def from_series_or_scalar(se):
+    try:
+        return DataArray.from_series(se)
+    except AttributeError:  # scalar case
+        return DataArray(se)
+
+
+def series_reduce(da, func, dim, **kwargs):
+    """ convert DataArray to pd.Series, apply pd.func, then convert back to
+    a DataArray. Multiple dims cannot be specified."""
+    if dim is None or da.ndim == 1:
+        se = da.to_series()
+        return from_series_or_scalar(getattr(se, func)(**kwargs))
+    else:
+        da1 = []
+        dims = list(da.dims)
+        dims.remove(dim)
+        d = dims[0]
+        for i in range(len(da[d])):
+            da1.append(series_reduce(da.isel(**{d: i}), func, dim, **kwargs))
+
+        if d in da.coords:
+            return concat(da1, dim=da[d])
+        return concat(da1, dim=d)
+
+
+@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'])
+@pytest.mark.parametrize('skipna', [False, True])
+@pytest.mark.parametrize('aggdim', [None, 'x'])
+def test_reduce(dim_num, dtype, dask, func, skipna, aggdim):
+
+    if aggdim == 'y' and dim_num < 2:
+        return
+
+    if dtype == np.bool_ and func == 'mean':
+        return  # numpy does not support this
+
+    if dask and not has_dask:
+        return
+
+    rtol = 1e-04 if dtype == np.float32 else 1e-05
+
+    da = construct_dataarray(dim_num, dtype, contains_nan=True, dask=dask)
+    axis = None if aggdim is None else da.get_axis_num(aggdim)
+
+    if dask and not skipna and func in ['var', 'std'] and dtype == np.bool_:
+        # TODO this might be dask's bug
+        return
+
+    if (LooseVersion(np.__version__) >= LooseVersion('1.13.0') and
+            da.dtype.kind == 'O' and skipna):
+        # Numpy < 1.13 does not handle object-type array.
+        try:
+            if skipna:
+                expected = getattr(np, 'nan{}'.format(func))(da.values,
+                                                             axis=axis)
+            else:
+                expected = getattr(np, func)(da.values, axis=axis)
+
+            actual = getattr(da, func)(skipna=skipna, dim=aggdim)
+            assert np.allclose(actual.values, np.array(expected), rtol=1.0e-4,
+                               equal_nan=True)
+        except (TypeError, AttributeError, ZeroDivisionError):
+            # TODO currently, numpy does not support some methods such as
+            # nanmean for object dtype
+            pass
+
+    # make sure the compatiblility with pandas' results.
+    actual = getattr(da, func)(skipna=skipna, dim=aggdim)
+    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)
+        expected = series_reduce(da, func, skipna=skipna, dim=aggdim, ddof=5)
+        assert_allclose(actual, expected, rtol=rtol)
+    else:
+        expected = series_reduce(da, func, skipna=skipna, dim=aggdim)
+        assert_allclose(actual, expected, rtol=rtol)
+
+    # make sure the dtype argument
+    if func not in ['max', 'min']:
+        actual = getattr(da, func)(skipna=skipna, dim=aggdim, dtype=float)
+        assert actual.dtype == float
+
+    # without nan
+    da = construct_dataarray(dim_num, dtype, contains_nan=False, dask=dask)
+    actual = getattr(da, func)(skipna=skipna)
+    expected = getattr(np, 'nan{}'.format(func))(da.values)
+    if actual.dtype == object:
+        assert actual.values == np.array(expected)
+    else:
+        assert np.allclose(actual.values, np.array(expected), rtol=rtol)
+
+
+@pytest.mark.parametrize('dim_num', [1, 2])
+@pytest.mark.parametrize('dtype', [float, int, np.float32, np.bool_, str])
+@pytest.mark.parametrize('contains_nan', [True, False])
+@pytest.mark.parametrize('dask', [False, True])
+@pytest.mark.parametrize('func', ['min', 'max'])
+@pytest.mark.parametrize('skipna', [False, True])
+@pytest.mark.parametrize('aggdim', ['x', 'y'])
+def test_argmin_max(dim_num, dtype, contains_nan, dask, func, skipna, aggdim):
+    # pandas-dev/pandas#16830, we do not check consistency with pandas but
+    # just make sure da[da.argmin()] == da.min()
+
+    if aggdim == 'y' and dim_num < 2:
+        return
+
+    if dask and not has_dask:
+        return
+
+    if contains_nan:
+        if not skipna:
+            # numpy's argmin (not nanargmin) does not handle object-dtype
+            return
+        if skipna and np.dtype(dtype).kind in 'iufc':
+            # numpy's nanargmin raises ValueError for all nan axis
+            return
+
+    da = construct_dataarray(dim_num, dtype, contains_nan=contains_nan,
+                             dask=dask)
+
+    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)
+    assert_allclose(actual.drop(actual.coords), expected.drop(expected.coords))
+
+
+def test_argmin_max_error():
+    da = construct_dataarray(2, np.bool_, contains_nan=True, dask=False)
+    with pytest.raises(ValueError):
+        da.argmin(dim='y')
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index 5a89627a..d018df8a 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -468,10 +468,17 @@ class VariableSubclassTestCases(object):
         a = self.cls('x', ['0', '1', '2'])
         b = self.cls('x', ['3', '4'])
         actual = Variable.concat([a, b], dim='x')
-        expected = Variable('x', np.arange(5).astype(str).astype(object))
+        expected = Variable('x', np.arange(5).astype(str))
         assert_identical(expected, actual)
-        assert expected.dtype == object
-        assert type(expected.values[0]) == str
+        assert actual.dtype.kind == expected.dtype.kind
+
+    def test_concat_mixed_dtypes(self):
+        a = self.cls('x', [0, 1])
+        b = self.cls('x', ['two'])
+        actual = Variable.concat([a, b], dim='x')
+        expected = Variable('x', np.array([0, 1, 'two'], dtype=object))
+        assert_identical(expected, actual)
+        assert actual.dtype == object
 
     def test_copy(self):
         v = self.cls('x', 0.5 * np.arange(10), {'foo': 'bar'})