New infer_intervals keyword for pcolormesh #1079

fmaussion · 2016-11-03T22:35:29Z

Addresses #781

@jhamman what do you think? I'm not sure if infer_interval_breaks is the best name for the keyword but I didn't come up with anything better than that. (maybe infer_coord_intervals?)

shoyer · 2016-11-04T01:01:23Z

I would call this just infer_intervals. breaks in the helper function name is referring to an implementation detail of how we store the intervals (by using breaks).

ocefpaf · 2016-11-04T09:12:47Z

Is infer_intervals creating coordinate bounds for plotting? It looks like it. What about when we do have cell boundaries in the dataset? (See http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#cell-boundaries). Will those be ignored and the inferred one used?

(Sorry if I am only making noise here and this does not make sense. In that case just ignore my comment.)

fmaussion · 2016-11-04T09:22:01Z

Is infer_intervals creating coordinate bounds for plotting? It looks like it.

Yes, it is a special case for pcolormesh though

What about when we do have cell boundaries in the dataset?

Currently the entire xarray workflow is built upon the implicit assumption that the coordinates are always at the grid-point center ( @shoyer correct me if I'm wrong). I'm not sure if the additional complexity added by cell boundaries is on the xarray devs priority list...

ocefpaf · 2016-11-04T09:33:42Z

Yes, it is a special case for pcolormesh though

A special case that I love. In the old Matlab day we had horrible hacks to get pcolor to work properly 😄

I'm not sure if the additional complexity added by cell boundaries is on the xarray devs priority list...

That is fine. Maybe a warning in the docs would be nice. So people can revert to a "manual" plotting to get the boundaries right.

fmaussion · 2016-11-04T10:03:04Z

In the old Matlab day we had horrible hacks to get pcolor to work properly

pcolormesh and imshow are indeed very consistent in xarray:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

da = xr.DataArray(np.arange(20).reshape(4, 5), dims=['lat', 'lon'],
                  coords = {'lat': np.linspace(0, 30, 4),
                            'lon': np.linspace(-20, 20, 5)})

f = plt.figure(figsize=(5, 8))

ax = plt.subplot(3, 1, 1, projection=ccrs.PlateCarree())
da.plot.imshow(ax=ax, transform=ccrs.PlateCarree())
lon, lat = np.meshgrid(da.lon.values, da.lat.values)
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_title('imshow')

ax = plt.subplot(3, 1, 2, projection=ccrs.PlateCarree())
da.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_title('pcolormesh')

ax = plt.subplot(3, 1, 3, projection=ccrs.PlateCarree())
da.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), infer_intervals=False)
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_title('pcolormesh - no intervals')

In most cases setting infer_intervals=False won't be a good idea, but maybe for some exotic projections (at the pole in @jhamman example) it might be easier this way (albeit a bit "wrong").

shoyer · 2016-11-04T15:08:34Z

We definitely ignore cell boundaries -- they don't (yet) have any place in the xarray data model.

It would be nice to have a built in "interval" data type, along with an IntervalIndex, that we could use for storing and doing lookups along interval axes. At one point I was working on this for pandas (see pandas-dev/pandas#7640), but then abandoned my effort when it seemed like it would be too much work to get it finished.

ocefpaf · 2016-11-04T15:26:25Z

We definitely ignore cell boundaries -- they don't (yet) have any place in the xarray data model.

Even though I would love to have that functionality I do not believe it is high priority. cell boundaries is a corner case, at least in my field of work, and 99% of the time it is OK to infer the intervals. Also, the name of the key word is quite clear: infer_intervals.

Maybe you could only add a note in the docs mentioning that the cell boundaries might exist?

jhamman · 2016-11-04T15:34:41Z

xarray/plot/plot.py

    """
    Pseudocolor plot of 2d DataArray

    Wraps matplotlib.pyplot.pcolormesh
    """

-    if not hasattr(ax, 'projection'):
+    if infer_intervals:


I have mixed feelings about this change. I'm not sure I like infer_intervals to default to True. This is going to break a bunch of our plotting code. My preference, and I'm open to other ideas at this point, would be to default to infer_intervals=None with something like this logic:

if (infer_intervals is None) and (hasattr(ax, 'projection')) infer_intervals = False warnings.warn("infer_intervals=None is deprecated, please specify either True or False", DeprecationWarning) if infer_intervals: ...

I'm not sure to understand... infer_intervals=True was the de-facto default until V0.7.2, and from the example above it is clear that setting it to False is almost never what you want to do. Can you provide an example where this is a good thing?

I agree your concerns about breaking old code though. How about keeping the current behavior for one more version, but issuing a warning that in the next version (probably 0.9.1 or 0.10), it is going to be set to True per default?

For my purposes this is going to imply that I'm going to use imshow() and not pcolormesh anymore for the time being.

Setting infer_intervals=False is exactly what you want if you have irregularly spaced coordinates (see example here: #781 (comment)).

I think something like the logic I showed above may be a good compromise. I'm happy to transition to always specifying infer_intervals=False we should be more transparent about the change (hence the warning along with a None for a default value. Would you be happier with this logic?

if (infer_intervals is None): if (hasattr(ax, 'projection')) infer_intervals = False warnings.warn("infer_intervals=None is deprecated, please specify either True or False", DeprecationWarning) else: infer_intervals = True if infer_intervals: ...

OK, with irregular coordinates none of the solutions is satisfying. infer_intervals=True crops one column, while infer_intervals=False crop one column and one line of data:

import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import xarray as xr lon, lat = np.linspace(-20, 20, 5), np.linspace(0, 30, 4) lon, lat = np.meshgrid(lon, lat) da = xr.DataArray(np.arange(20).reshape(4, 5), dims=['y', 'x'], coords = {'lat': (('y', 'x'), lat), 'lon': (('y', 'x'), lon)}) f = plt.figure() ax = plt.subplot(2, 1, 1, projection=ccrs.PlateCarree()) da.plot.pcolormesh('lon', 'lat', ax=ax, transform=ccrs.PlateCarree(), infer_intervals=True) ax.scatter(lon, lat, transform=ccrs.PlateCarree()) ax.coastlines() ax.set_title('pcolormesh') ax = plt.subplot(2, 1, 2, projection=ccrs.PlateCarree()) da.plot.pcolormesh('lon', 'lat', ax=ax, transform=ccrs.PlateCarree(), infer_intervals=False) ax.scatter(lon, lat, transform=ccrs.PlateCarree()) ax.coastlines() ax.set_title('pcolormesh - no intervals') plt.show()

I'll try to have a closer look later next week

fmaussion · 2016-11-06T17:22:14Z

After giving this a few more thoughts:

if the coordinates are 1D (i.e. the grid is regular in the scipy sense: the coordinates don't have to be evenly spaced), infer_intervals should be set to True by default. This should make everybody happy (including @jhamman and me), and this will cover a huge majority of the geoscientific grids.
if the coordinates are 2D, what xarray was doing until now was wrong. _infer_interval_breaks was inferring intervals on the first axis only and ignoring the second, so regardless of whether we were plotting on a projection or not, the plot output was cropped (see example above)

So what do we do with this? @jhamman 's suggestion is to set infer_intervals to False, which is equivalent to saying: "this plot is going to be a bit wrong but this is not xarray's problem". I guess that it's fine to do so, but we'd have to be happy with the following:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

lon, lat = np.linspace(-20, 20, 5), np.linspace(0, 30, 4)
lon, lat = np.meshgrid(lon, lat)
lon += lat/10
lat += lon/10
da = xr.DataArray(np.arange(20).reshape(4, 5), dims=['y', 'x'],
                  coords = {'lat': (('y', 'x'), lat),
                            'lon': (('y', 'x'), lon)})

ax = plt.subplot(1, 2, 1, projection=ccrs.PlateCarree())
da.plot.pcolormesh('lon', 'lat', ax=ax, transform=ccrs.PlateCarree(), infer_intervals=False)
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_title('2d irreg case - pcolor')

ax = plt.subplot(1, 2, 2, projection=ccrs.PlateCarree())
da.plot.contourf('lon', 'lat', ax=ax, transform=ccrs.PlateCarree())
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_title('2d irreg case - contourf')

Note that one line and one column of data disappeared (so that the colorbar range doesn't correspond to what we see) and the coordinates are not understood as being in the cell center. The advantage of this choice are twofold:

it's kind of coherent with how contourf would interpret the coordinates (see example above)
if we agree on "1D = infer per default" and "2D = don't infer per default", we might even get rid of the new keyword altogether.

Let's call this option "Option A".

[edited]

I've been playing around with an "Option B", which includes making _infer_interval_breaks accept an axis kwarg. Regardless of the implementation, we could reach the following result:

[\edited]

This plot is closer to what I think the grid actually represents, but I am biased towards my usage of xarray (gridded geoscientific models). The disadvantages of option B are:

it adds a bit of complexity
if you think about it, the problem of plotting irregular grids is not really an xarray problem. If there was a perfect way to handle an irregular grids without having bounded coordinates, I guess that matplotlib would have a tool for it.

I currently tend to go for option A, since option B is a bit awkward and users who have bounded coordinated available somewhere (as suggested by @ocefpaf ) will have to implement their own plotting routine anyway.

fmaussion · 2016-11-07T13:55:56Z

I've come up with a cleaner implementation which doesn't require scipy and produces a very accurate plot:

I'm going to clean all this later today and update my PR, I think that we can come to something nice.

@jhamman : I wonder if the problem you are having in your example is related to flaoting point issues at the pole. Would you mind sending me the file and the code that produces the plot so that I can have a look at it?

jhamman · 2016-11-07T14:23:35Z

You can get data on the same grid (different variable, same coords) from the sample data ds = xr.tutorial.load_dataset('rasm').isel(time=0). From there, my example should work.

fmaussion · 2016-11-07T18:33:20Z

@jhamman @shoyer this is ready for another review (the test failure seems unrelated).

@ocefpaf I added a note to the docs

@jhamman OK now I get why it clearly doesn't make sense to infer anything in your dataset (plot below). Just out of curiosity: does RASM have a map projection on which it has regular coordinates, too?

fmaussion · 2016-11-07T18:43:24Z

Actually, you wouldn't like to infer intervals on these coordinates at all, regardless if you are on a map or not. Should I remove the if hasattr(ax, 'projection') altogether?

fmaussion · 2016-11-08T18:15:45Z

Should I remove the if hasattr(ax, 'projection') altogether?

Finally I think it's ok to keep it as it is now, since the data above is likely to always be plotted on a map anyways.

rabernat · 2016-11-08T18:31:16Z

This looks like a really nice PR.

There is some intersection with my "multidimensional coordinates" tutorial:
http://xarray.pydata.org/en/stable/examples/multidimensional-coords.html

At the least, we should standardize the terminology in the docs. Are these "irregular coordinates" or "multidimensional coordinates"? I don't particularly care, but we should be consistent.

fmaussion · 2016-11-08T18:36:43Z

@rabernat thanks for the link! I'll update the docs with a link to your example and with "multidimensional" instead of "irregular".

fmaussion · 2016-11-08T19:39:40Z

@rabernat : couldn't your tutorial also use the ipython directive and run at build time instead of using the parsed literals? I see that the other examples are also static, is there a reason for this?

rabernat · 2016-11-08T19:48:01Z

is there a reason for this?

My laziness / ignorance of the documentation build process.

fmaussion · 2016-11-08T19:51:46Z

OK. Let's see what the gurus say about this, maybe there is a reason to keep it like this

shoyer · 2016-11-08T20:16:58Z

Yes, could certainly run this at doc build time instead of using a notebook. The advantage of the notebook is that it's easier to edit/assemble, and also easier to download and adjust.

In an ideal world, I think we would check in a notebook with cleared output and then run it and convert to HTML at doc build time. Looks like that's exactly what nbsphinx does if anyone wants to take a go at setting that up. Could be a much nicer way to write the docs than RST.

fmaussion · 2016-11-08T20:25:21Z

Note that in that case, multidimensional-coords and monthly-means are already formatted in rst.

shoyer · 2016-11-08T20:26:08Z

Right, we used a static ipynb -> rst converter

jhamman

I'm fine with this fix. I know my suggestions lead to a bunch of complicated logic but I think this is the best solution for now. Thanks for putting in the hard work to get here.

shoyer · 2016-11-10T02:23:16Z

xarray/plot/plot.py

-    last = coord[-1] + deltas[-1]
-    return np.r_[[first], coord[:-1] + deltas, [last]]
+    if axis == 0:
+        deltas = 0.5 * (coord[1:, ...] - coord[:-1, ...])


You could write an axis agnostic version of this with np.diff, np.take and np.concatenate, e.g.,

deltas = 0.5 * np.diff(coord, axis) first = np.take(coord, [0], axis) last = np.take(coord, [-1], axis) trim_last = tuple(slice(None, -1) if n == axis else slice(None) for n in range(coord.ndim)) return np.concatenate([first, coord[trim_last] + deltas, last], axis)

Maybe a slight improvement? (Up to you on whether to use it)

shoyer · 2016-11-10T02:23:37Z

xarray/plot/plot.py

-    first = coord[0] - deltas[0]
-    last = coord[-1] + deltas[-1]
-    return np.r_[[first], coord[:-1] + deltas, [last]]
+    if axis == 0:


Add an assert or error to verify that axis is one of the expected values.

shoyer · 2016-11-10T02:24:38Z

xarray/plot/plot.py

+        else:
+            # we have to infer the intervals on both axes
+            x = _infer_interval_breaks(x, axis=1)
+            x = _infer_interval_breaks(x, axis=0)


fmaussion · 2016-11-10T16:21:28Z

I had to modify your code a little bit for this to work, but this is now accepting any axis. Quite a useful function actually!

ocefpaf · 2016-11-10T22:19:19Z

xarray/plot/plot.py

+    last = np.take(coord, [-1], axis=axis) + np.take(deltas, [-1], axis=axis)
+    trim_last = tuple(slice(None, -1) if n == axis else slice(None)
+                      for n in range(coord.ndim))
+    return np.concatenate([first, coord[trim_last] + deltas, last], axis=axis)


Nifty! Loved the use of axis here!

shoyer · 2016-11-10T22:27:29Z

xarray/test/test_plot.py

@@ -112,6 +112,17 @@ def test__infer_interval_breaks(self):
        self.assertArrayEqual(pd.date_range('20000101', periods=4) - np.timedelta64(12, 'h'),
                              _infer_interval_breaks(pd.date_range('20000101', periods=3)))

+        # make a bounded 2D array that we will center and re-infer
+        xref, yref = np.meshgrid(np.arange(6), np.arange(5))
+        cx = (xref[1:, 1:] + xref[:-1, :-1]) / 2


On Python 2, I think this will still have integer dtype. Can you divide by 2.0 instead?

I added a from __future__ import division instead

shoyer · 2016-11-10T22:39:40Z

xarray/test/test_plot.py

@@ -1,3 +1,5 @@
+from __future__ import division


Can you add this to the top of literally every single Python file in xarray?

That will keep this easier to reason about.

Actually:

from __future__ import absolute_import from __future__ import division from __future__ import print_function

Yes OK. shouldn't I do this in a separate PR? I'll do it tomorrow though, it's getting late in Innsbruck...

shoyer · 2016-11-10T22:54:15Z

Yes, thanks!

On Thu, Nov 10, 2016 at 2:43 PM, Fabien Maussion [email protected]
wrote:

@fmaussion commented on this pull request.

In xarray/test/test_plot.py #1079:

@@ -1,3 +1,5 @@
+from future import division

Yes OK. shouldn't I do this in a separate PR? I'll do it tomorrow though,
it's getting late in Innsbruck...

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#1079, or mute the thread
https://github.com/notifications/unsubscribe-auth/ABKS1vqdXfCSLnXeKwewm2EH6rqzvp1fks5q854MgaJpZM4Ko_0w
.

fmaussion added 2 commits November 3, 2016 23:33

fixes pydata#781

8584b91

typo

182280a

rename keyword

f8953fd

the kwargs should also be passed over now

26dc736

jhamman reviewed Nov 4, 2016

View reviewed changes

fmaussion changed the title ~~New infer_interval_breaks keyword for pcolormesh~~ New infer_intervals keyword for pcolormesh Nov 4, 2016

rabernat mentioned this pull request Nov 4, 2016

N + 1 sized grids for X and Y MITgcm/xmitgcm#15

Open

infer 2d coords, new chapter in docs

a40ab58

revert to original cartopy test

923af9d

update docs

95513dc

jhamman approved these changes Nov 9, 2016

View reviewed changes

shoyer reviewed Nov 10, 2016

View reviewed changes

other version of _infer_intervals_breaks

d4ff6f9

ocefpaf reviewed Nov 10, 2016

View reviewed changes

shoyer reviewed Nov 10, 2016

View reviewed changes

py2 division

b9ec680

shoyer reviewed Nov 10, 2016

View reviewed changes

shoyer merged commit 8ded250 into pydata:master Nov 10, 2016

fmaussion mentioned this pull request Nov 11, 2016

py2 compat header #1108

Merged

fmaussion mentioned this pull request Jul 11, 2017

Allow DataArray to hold cell boundaries as coordinate variables #1475

Open

fmaussion mentioned this pull request Jan 24, 2018

bug: 2D pcolormesh plots are wrong when coordinate is not ascending order #1852

Closed

New infer_intervals keyword for pcolormesh #1079

New infer_intervals keyword for pcolormesh #1079

Conversation

fmaussion commented Nov 3, 2016

shoyer commented Nov 4, 2016

ocefpaf commented Nov 4, 2016

fmaussion commented Nov 4, 2016

ocefpaf commented Nov 4, 2016

fmaussion commented Nov 4, 2016

shoyer commented Nov 4, 2016

ocefpaf commented Nov 4, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fmaussion commented Nov 6, 2016 • edited Loading

fmaussion commented Nov 7, 2016

jhamman commented Nov 7, 2016

fmaussion commented Nov 7, 2016

fmaussion commented Nov 7, 2016

fmaussion commented Nov 8, 2016

rabernat commented Nov 8, 2016

fmaussion commented Nov 8, 2016

fmaussion commented Nov 8, 2016 • edited Loading

rabernat commented Nov 8, 2016

fmaussion commented Nov 8, 2016

shoyer commented Nov 8, 2016

fmaussion commented Nov 8, 2016 • edited Loading

shoyer commented Nov 8, 2016

jhamman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fmaussion commented Nov 10, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shoyer commented Nov 10, 2016

@fmaussion commented on this pull request.

fmaussion commented Nov 6, 2016 •

edited

Loading

fmaussion commented Nov 8, 2016 •

edited

Loading

fmaussion commented Nov 8, 2016 •

edited

Loading