Repository: pydata/xarray

Issue:
bug: 2D pcolormesh plots are wrong when coordinate is not ascending order
#### Code Sample, a copy-pastable example if possible

```python
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr


x = np.arange(10)
y = np.arange(20)

np.random.shuffle(x)

x = xr.DataArray(x, dims=['x'], coords={'x': x})
y = xr.DataArray(y, dims=['y'], coords={'y': y})

z = x + y
z_sorted = z.isel(x=np.argsort(x.values))


# make plot
fig, axs= plt.subplots(1, 2, figsize=(6,3))
z_sorted.plot(ax=axs[0])
axs[0].set_title("X is sorted")

z.plot(ax=axs[1])
axs[1].set_title("X is not unsorted")
plt.tight_layout()
```
#### Problem description

Sometime the coordinates in an xarray dataset are not always sorted in ascending order. I recently had an issue where the time coordinate of a 2D datasets was scrambled, so calling `x.plot` gave very strange results. In my opinion, `x.plot` should probably sort the data along the coordinates, or at least provide a warning if the coordinates are unsorted.

#### Expected Output
Here is the image generated by the snippet above:

![image](https://user-images.githubusercontent.com/1386642/35318676-36d136a0-0091-11e8-9c94-bce67f4d5991.png)

The left and right panels should be the same.

# Paste the output here xr.show_versions() here

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 3.6.2.final.0
python-bits: 64
OS: Darwin
OS-release: 16.0.0
machine: x86_64
processor: i386
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8

xarray: 0.10.0+dev50.ga988dc2
pandas: 0.20.3
numpy: 1.13.1
scipy: 0.19.1
netCDF4: 1.3.1
h5netcdf: 0.5.0
Nio: None
zarr: None
bottleneck: 1.2.1
cyordereddict: None
dask: 0.15.2
distributed: 1.18.3
matplotlib: 2.0.2
cartopy: None
seaborn: 0.8.0
setuptools: 36.5.0.post20170921
pip: 9.0.1
conda: 4.3.29
pytest: 3.2.1
IPython: 6.1.0
sphinx: 1.6.3
</details>


Discussion:
fmaussion: Thanks a lot @nbren12 ! Indeed `_infer_interval_breaks` is probably the right place to add this warning.

A few comments:
- the tests are now failing, likely because of a type issue with time coordinates
- coordinates don't have to be ascending for pcolormesh to work (see example below). Monotonically ascending and descending coordinates should work.
- we need a test for this new behavior ([asserting warnings](https://docs.pytest.org/en/3.0.2/recwarn.html))
- we follow the pep8 conventions, and I think your lines are longer than the 79 character rule

```python
import numpy as np
import xarray as xr
x = np.arange(5)[::-1]
y = np.arange(6)[::-1]
x = xr.DataArray(x, dims=['x'], coords={'x': x})
y = xr.DataArray(y, dims=['y'], coords={'y': y})
z = x + y
z.plot();
```

![index](https://user-images.githubusercontent.com/10050469/35766253-f1574730-08d4-11e8-9def-7a2f1fa13965.png)

nbren12: Thanks for your help @fmaussion. Hopefully the commit I just pushed fixes the failing tests. I will work on adding a test and fixing the formatting.
nbren12: Ok. I think everything is ready now.
nbren12: @shoyer That would work with me. Is there any chance people would want to
make heatmaps involving categorical variables though?

If we do decide to raise an error, why not go one step further and just
sort the coordinates automatically?

On Sun, Feb 4, 2018 at 4:00 PM, Stephan Hoyer <notifications@github.com>
wrote:

> *@shoyer* commented on this pull request.
> ------------------------------
>
> In xarray/plot/plot.py
> <https://github.com/pydata/xarray/pull/1885#discussion_r165863769>:
>
> > @@ -750,6 +767,13 @@ def _infer_interval_breaks(coord, axis=0):
>             [ 2.5,  3.5,  4.5]])
>      """
>      coord = np.asarray(coord)
> +
> +    if not _is_monotonic(coord, axis=axis):
> +        warnings.warn("The input coordinate is not sorted in increasing order "
>
> Rather an a warning, why not make this an error? I don't see any use-cases
> for 2d plots with non-monotonic coordinates. With the current version of
> xarray, these plots always end up wrong in some way, either by not plotting
> everything or with bad axis labels.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/pydata/xarray/pull/1885#pullrequestreview-93863697>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/ABUokuMwNBc8JRgafrF7HNmGfgzxi64Qks5tRkSbgaJpZM4R4EWU>
> .
>

shoyer: > Is there any chance people would want to make heatmaps involving categorical variables though?

Certainly, but the currently version of 2D plots in xarray doesn't work for these!

If they want to make plots involving categorical variables, I would currently recommend using `seaborn.heatmap` (we can mention it in the error message, too).

> If we do decide to raise an error, why not go one step further and just
sort the coordinates automatically?

One reason is that it's not obvious if they would like increasing or decreasing coordinates.
nbren12: >  One reason is that it's not obvious if they would like increasing or decreasing coordinates.

For me at least, `plt.pcolormesh` automatically displays in increasing order even if one of the input arrays is sorted in descending order. This happens all the time with meteorological data available in pressure coordinates (pressure goes down with height). I usually have to manually call plt.gca().invert_yaxis() to flip the y axis. Does `xarray.plot`'s behave differently?
nbren12: That’s a good point. Probably not.
On Tue, Feb 6, 2018 at 5:25 AM Fabien Maussion <notifications@github.com>
wrote:

> *@fmaussion* commented on this pull request.
> ------------------------------
>
> In xarray/plot/plot.py
> <https://github.com/pydata/xarray/pull/1885#discussion_r166295479>:
>
> > @@ -741,6 +741,23 @@ def contourf(x, y, z, ax, **kwargs):
>      return primitive
>
>
> +def _is_monotonic(coord, axis=0):
> +    """
> +    >>> _is_monotonic(np.array([0, 1, 2]))
> +    True
> +    >>> _is_monotonic(np.array([2, 1, 0]))
> +    True
> +    >>> _is_monotonic(np.array([0, 2, 1]))
> +    False
> +    """
> +    if coord.shape[axis] < 3:
> +        return True
> +    else:
> +        delta = np.diff(coord, axis=axis)
> +        # This ensures all the deltas have the same sign
> +        return np.abs(delta).sum() == np.abs(delta.sum())
>
> Is this robust wrt floating point accuracy?
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/pydata/xarray/pull/1885#pullrequestreview-94340385>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/ABUokjAvEYLAcFZ78V0_Ub_RYkZfDKuAks5tSFLlgaJpZM4R4EWU>
> .
>

fmaussion: Which, I guess, illustrates our point quite well: while it's easy for users to find out when their data is messy, it's much harder to find a robust and logical way to handle all cases for a general library like xarray. I guess it's the reason why matplotlib also doesn't warn or error when unordered data is given to pcolormesh.

Don't get me wrong: I think it would be valuable to be able to provide a better "service" than matplotlib. Probably we should even error earlier, at the 2d plot level (and not only at the pcolormesh level like it is now). It's just going to be a bit hard because of all the special cases (dtypes, 2d coords, and map projections moslty).

Ordering "per default" might be the most straightforward way, but:
- what about performance issues?
- I don't like the idea of silently doing this with user data

I'll let @shoyer decide on what's best here. I am slightly in favor of doing nothing and giving the responsibility to the user. 


nbren12: There. I think this code should work for all dtypes and 2D coords. 
fmaussion: Thanks! This looks good.
nbren12: Yah. I knew my solution before was probably too cute.

On Tue, Feb 6, 2018 at 11:14 AM, Fabien Maussion <notifications@github.com>
wrote:

> Thanks! This looks good.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/pydata/xarray/pull/1885#issuecomment-363532796>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/ABUokvVDiWCnX_JPa5TsDtW8JtQxdTDEks5tSKSZgaJpZM4R4EWU>
> .
>

fmaussion: Sorry @nbren12 there is still a pep8 error somewhere, and you can freely make it an error as @shoyer suggests.

Note that currently the test for monotonic coords  is for pcolormesh only, while it could be for all 2d plots as well. What I like with the current implementation is that the check happens in `_infer_interval_breaks`, which is only called when certain conditions are met (1D coords, or 2D coords without transform kwarg). If this check lives in the common 2d function the same conditions should be applied.

If we start to check input for sanity we might as well raise an error when coordinates are not regularly spaced in the imshow case, although I can imagine this might be more difficult to check.

Thoughts?
nbren12: > Sorry @nbren12 there is still a pep8 error somewhere, and you can freely make it an error as @shoyer suggests.

Hopefully, the commit I just pushed fixes this. I had some global flake8 settings that were messing with  my local linting.

> Note that currently the test for monotonic coords is for `pcolormesh` only, while it could be for all 2d plots as well.

I don't think I was having any issues with `contour` or `contourf` when I first opened the issue, but it probably does break `imshow`. 

> If we start to check input for sanity we might as well raise an error when coordinates are not regularly spaced in the imshow case.

Maybe we could leave this to a later PR.
fmaussion: Sorry, our new flake8 tool can be quite strict ;)
nbren12: no problem. I have always preferred putting operators on new lines, though. Didn't realize that was against pep8.  oh well 🤷‍♂️ . 
nbren12: I just rebased onto master.

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index fe8332c3..6018a184 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -106,6 +106,8 @@ Enhancements
 
 Bug fixes
 ~~~~~~~~~
+- Support indexing with a 0d-np.ndarray (:issue:`1921`).
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
 - Added warning in api.py of a netCDF4 bug that occurs when
   the filepath has 88 characters (:issue:`1745`).
   By `Liam Brannigan <https://github.com/braaannigan>`_.
diff --git a/xarray/core/variable.py b/xarray/core/variable.py
index 14fc3480..267dc02c 100644
--- a/xarray/core/variable.py
+++ b/xarray/core/variable.py
@@ -463,10 +463,14 @@ class Variable(common.AbstractArray, utils.NdimSizeLenMixin):
         key = self._item_key_to_tuple(key)  # key is a tuple
         # key is a tuple of full size
         key = indexing.expanded_indexer(key, self.ndim)
-        # Convert a scalar Variable as an integer
+        # Convert a scalar Variable to an integer
         key = tuple(
             k.data.item() if isinstance(k, Variable) and k.ndim == 0 else k
             for k in key)
+        # Convert a 0d-array to an integer
+        key = tuple(
+            k.item() if isinstance(k, np.ndarray) and k.ndim == 0 else k
+            for k in key)
 
         if all(isinstance(k, BASIC_INDEXING_TYPES) for k in key):
             return self._broadcast_indexes_basic(key)
diff --git a/xarray/plot/plot.py b/xarray/plot/plot.py
index 97fee5f0..162a0c23 100644
--- a/xarray/plot/plot.py
+++ b/xarray/plot/plot.py
@@ -745,6 +745,26 @@ def contourf(x, y, z, ax, **kwargs):
     return primitive
 
 
+def _is_monotonic(coord, axis=0):
+    """
+    >>> _is_monotonic(np.array([0, 1, 2]))
+    True
+    >>> _is_monotonic(np.array([2, 1, 0]))
+    True
+    >>> _is_monotonic(np.array([0, 2, 1]))
+    False
+    """
+    if coord.shape[axis] < 3:
+        return True
+    else:
+        n = coord.shape[axis]
+        delta_pos = (coord.take(np.arange(1, n), axis=axis) >=
+                     coord.take(np.arange(0, n-1), axis=axis))
+        delta_neg = (coord.take(np.arange(1, n), axis=axis) <=
+                     coord.take(np.arange(0, n-1), axis=axis))
+        return np.all(delta_pos) or np.all(delta_neg)
+
+
 def _infer_interval_breaks(coord, axis=0):
     """
     >>> _infer_interval_breaks(np.arange(5))
@@ -754,6 +774,15 @@ def _infer_interval_breaks(coord, axis=0):
            [ 2.5,  3.5,  4.5]])
     """
     coord = np.asarray(coord)
+
+    if not _is_monotonic(coord, axis=axis):
+        raise ValueError("The input coordinate is not sorted in increasing "
+                         "order along axis %d. This can lead to unexpected "
+                         "results. Consider calling the `sortby` method on "
+                         "the input DataArray. To plot data with categorical "
+                         "axes, consider using the `heatmap` function from "
+                         "the `seaborn` statistical plotting library." % axis)
+
     deltas = 0.5 * np.diff(coord, axis=axis)
     if deltas.size == 0:
         deltas = np.array(0.0)


Test patch:
diff --git a/xarray/tests/test_plot.py b/xarray/tests/test_plot.py
index 479b1dce..46410cd5 100644
--- a/xarray/tests/test_plot.py
+++ b/xarray/tests/test_plot.py
@@ -178,6 +178,10 @@ class TestPlot(PlotTestCase):
         np.testing.assert_allclose(xref, x)
         np.testing.assert_allclose(yref, y)
 
+        # test that warning is raised for non-monotonic inputs
+        with pytest.raises(ValueError):
+            _infer_interval_breaks(np.array([0, 2, 1]))
+
     def test_datetime_dimension(self):
         nrow = 3
         ncol = 4
diff --git a/xarray/tests/test_variable.py b/xarray/tests/test_variable.py
index d018df8a..c8f74762 100644
--- a/xarray/tests/test_variable.py
+++ b/xarray/tests/test_variable.py
@@ -627,6 +627,12 @@ class VariableSubclassTestCases(object):
         v_new = v[np.array([0])[0]]
         assert_array_equal(v_new, v_data[0])
 
+        v_new = v[np.array(0)]
+        assert_array_equal(v_new, v_data[0])
+
+        v_new = v[Variable((), np.array(0))]
+        assert_array_equal(v_new, v_data[0])
+
     def test_getitem_fancy(self):
         v = self.cls(['x', 'y'], [[0, 1, 2], [3, 4, 5]])
         v_data = v.compute().data