new shading default can mess up plots #1638
Comments
|
Yes, Matplotlib is deprecating flat shading because it was silently dropping the final row/column of data when the shapes were equal. With the update, See here for a good description: Unfortunately, this dataset is not strictly increasing and so the difference between coordinates jumps from 360 to 0 in the middle of the xc array: You can take Cartopy away completely and just plot in x/y coordinates and the problem is evident there too. |
|
Then this is an upstream issue - thanks. Now the quest begins to find out if there is already a mpl issue on this. Feel free to close. I add an example that works without xarray: import numpy as np
import matplotlib.pyplot as plt
lon = np.concatenate((np.arange(180, 360, 5), np.arange(0, 180, 5)))
lat = np.arange(90, 0, -1)
data = np.random.randn(*lat.shape + lon.shape) + lon
# set some data to nan to make the problem more visible
data[::5, 35] = np.NaN
plt.pcolormesh(lon, lat, data, shading="auto") |
|
Ahh but in pure matplotlib this does not work with rasm.isel(time=1).Tair.plot.pcolormesh.plot(x="xc", y="yc", shading="flat")
rasm.isel(time=1).Tair.plot.pcolormesh.plot(x="xc", y="yc", shading="nearest")So mpl can probably assume that the coordinates are monotonic. |
|
So, the coordinates don't have to be monotonic, but the cells will all be "connected" in a sense. If you open up an issue, feel free to ping me on it. I think that a warning should probably be added to Matplotlib for the new 'auto' case when things are not strictly increasing/decreasing in the coordinate arrays and they try to find midpoints. |
|
I closed the upstream issue just so we can hash this out here. Of course Matplotlib is willing to make changes, particularly if the new default is going to end up being problematic. I'm having trouble seeing what the way forward here is. The dataset above has a wrapping issue with no Cartopy projection in either "flat" or "nearest". Looking at the data set, that seems to be because the wrap occurs for different rows in different columns of I assume that flat works with NorthPolarStereo because the data set was developed for that projection. If we add a warning in Matplotlib, is it really just for "nearest" that there is an issue? I guess the "bad" thing we do for nearest is add a data point, which in this case is incorrect. I'm wondering if we should really do the grid re-centring in axes space and then transform back to data space... |
|
The thing I am still confused by is how the original xarray example ever worked in the previous versions because the wrap point was in the middle of the array and not one the values that was clipped. I'm curious if Cartopy previously did just fine transforming the coordinates and then identified those cells as wrapped and masked them out. So, Cartopy was essentially hiding the previous pcolormesh behavior. The original example in pure MPL still is not correct even with MPL 3.2 import matplotlib.pyplot as plt
import xarray as xr
rasm = xr.tutorial.load_dataset("rasm")
ax = plt.subplot(111)
rasm.isel(time=1).Tair.plot.pcolormesh(
ax=ax, x="xc", y="yc", shading='flat')
plt.show() |
|
The grid in that example looks like this in Cartesian co-ordinates (every 10th row and column)
If you form a quadrilateral in cartesian co-ordinates (default no cartopy) you have the ones near the edge with at least one corner that wraps. Presumably when you plot these from `NorthPolarStereo' there is no wrap. |
|
Just as prior art, xarray has |
|
I gave this a bit of thought and am not sure where we should fix this, either the MPL or Cartopy side. Cartopy is currently using a private method from MPL I think the easiest way out is to add a keyword to dX = np.diff(X, axis=1)
if wrap is not None:
# Account for angle wraps
dX = (dX + wrap/2) % wrap - wrap/2
dX = dX/2.Any other ideas on how we can keep Cartopy and MPL in sync here? |
|
@greglucas There are two types of wraps here, if I understand correctly. The simplest is what happens at -179 to +179 where there is mathematically a wrap, but after applying the projection there is no wrap i.e. in screen space the two data points can be next to each other (depending on central longitude of course). I think matplotlib could handle that wrap by doing the grid re-interpolation (or recentering) in screen space (i.e. after the projection transform has been applied) rather than before, which it does now. This would subtly skew existing grids for non-linear projections, but I don't think it would be too bad. Your proposed The second type of wrapping is due to the user not giving you the data properly sorted into physical space. eg if the user gives you data from -180 to +180, but then makes the central longitude 80 degrees. The axes now go from -100 to -100, and data at -99 is not contiguous on screen with -101, but it is in the data set. In that case I think cartopy is on its own because matplotlib can't re-order the data or insert masks to create gaps. |
|
Correct, I think we are on the same page with multiple wraps. Cartopy handles the "screen"/"projection" wraps currently. It can handle the points [181, 360] and map those to the negatives [-179, 0] and masks the the values that cross the wrap boundary and internally redraws those quads. The part that is messing Cartopy up currently is handling inserting the new points to make the quad edges that aren't correct in the data-space. So, really all we need to do is add in proper data-space points and then Cartopy will handle the rest. That is the easy fix and will work with data going The harder case is actually if a user decides to put their data in this order: I would have to give some more thought to how we could identify and mask quads if we do the interpolation in screen space instead. I am generally more of a proponent of doing the interpolations in the space that the data is given, not the one it is being mapped to. If a user gives us linearly spaced points, but wants to display them on a log-scaled axis do we add the bin edges in log-space (screen coordinates) or linear-space (data coordinates)? I'm not sure there is a correct answer here, really it probably boils down to preference. |
I generally agree with you, but keep in mind its just a slightly different dx we are talking about - the quad centres will still be spaced as before, just their extents would be a bit different. To me that is better than completely breaking if a non-linear projection is used. If a user is fussy that the quad edges line up exactly 0.5 of the cell width in data space, then they can specify the quad edges explicitly and not rely on whatever we do. |
|
I just pushed up an option to handle this on Cartopy's side in #1646. That handles defining bin edges in data space before getting to Matplotlib. I think that would be agnostic of whether you change MPL to do screenspace interpolation or not. I guess that I find this way cleaner personally, but it does add some of the MPL grid interpolation code into Cartopy now. |
|
Without checking exactly what you did I think that's a great solution for downstream: don't count on matplotlib generic edge interpolation if you don't need to. |
|
OK, so for clarity Matplotlib has decided to not deal with fancy re-interpolation. If users encounter non-monotonic grid, indicative of a possible wrap issue, they now get a warning and are asked to specify the bin edges rather than the centers. (Thanks matplotlib/matplotlib#18398 @greglucas!). If we continue to get feedback that this is inadequate we may consider trying to do something fancier like interpolation in display space, but it wasn't clear if people want quadrilaterals to be square in display space of data space, so we decided to punt for matplotlib 3.3.2. |
|
Sounds good. I think the new warning message is pretty clear now. And hopefully it is only a small number of people that will ever hit that path in the first place... |
|
Hmm, I think a lot of people will hit it. I plot pcolormesh all the time that double back on themselves. e.g. we have an underwater glider trying to drive west, but its getting pushed around some by currents, so sometimes goes backwards. A pcolormesh is perfectly reasonable in this case, and its perfectly reasonable to not have the x co-ordinate be monotonic. |
|
Interesting use case, neat. I generally bin/grid my data so it always ends up in some sorted order. I guess you're just OK with overlapping quads and hiding the 'lower' one, assuming your most recent data is best? |
|
One hopes the data aren't that different, and of course you can always plot in time if you want to unfold it. Binning the data is OK for some applications, but involves interpolation and/or averaging which is not always a luxury the data can afford. |
Description
The new shading default
cartopy/lib/cartopy/mpl/geoaxes.py
Line 1669 in 08db54f
can lead to artifacts in the plot which do not happen with
shading="flat". This happens for coordinates that cross the zero line. As far as I can see matplotlib plans to entirely remove this option?Code to reproduce
I was too lazy to create an example without xarray as dependency - let me know if you need one:
Artifacts:
No artifacts:
Full environment definition
Operating system
linux
Cartopy version
0.18.0
conda list
pip list
The text was updated successfully, but these errors were encountered: