issue_comments

Yes, thanks!

On Thu, Nov 10, 2016 at 2:43 PM, Fabien Maussion notifications@github.com wrote:

@fmaussion commented on this pull request.

In xarray/test/test_plot.py https://github.com/pydata/xarray/pull/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 https://github.com/pydata/xarray/pull/1079, or mute the thread https://github.com/notifications/unsubscribe-auth/ABKS1vqdXfCSLnXeKwewm2EH6rqzvp1fks5q854MgaJpZM4Ko_0w .

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!

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

Right, we used a static ipynb -> rst converter

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.

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

is there a reason for this?

My laziness / ignorance of the documentation build process.

@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 thanks for the link! I'll update the docs with a link to your example and with "multidimensional" instead of "irregular".

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.

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.

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?

@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?

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.

After giving this a few more thoughts: 1. 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. 2. 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:

``` python 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.

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?

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?

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 https://github.com/pandas-dev/pandas/issues/7640), but then abandoned my effort when it seemed like it would be too much work to get it finished.

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

pcolormesh and imshow are indeed very consistent in xarray:

``` python 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").

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 :smile:

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.

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...

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.)

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).