pcolormesh: shading='nearest' and non-monotonic coordinates

@greglucas

Bug report

Bug summary

#16258 deprecates shading="flat" when data and coordinates have the same shape and introduces shading="nearest" which interpolates the coordinates. However, the logic to interpolate the coordinates breaks down if they are not monotonically increasing.

matplotlib/lib/matplotlib/axes/_axes.py

Line 5684 in 8511771

def _interp_grid(X):

Of course non-monotonic coordinates are a bit exotic but they make sense for geographic data.

Code for reproduction

import numpy as np

def _interp_grid(X):
    # helper for below
    if np.shape(X)[1] > 1:
        dX = np.diff(X, axis=1)/2.
        X = np.hstack((X[:, [0]] - dX[:, [0]],
                       X[:, :-1] + dX,
                       X[:, [-1]] + dX[:, [-1]]))
    return X

x = np.array([1, 2, 0])
y = np.array([1, 2])

X, Y = np.meshgrid(x, y)

_interp_grid(X)

Actual outcome

Results in

array([[ 0.5,  1.5,  1. , -1. ],
       [ 0.5,  1.5,  1. , -1. ]])

While having cooridinates [1, 2, 0] may seem just as wrong as coordinates [0.5, 1.5, 1, -1] the latter can cause some real-world problems in conjunction with cartopy (SciTools/cartopy#1638):

Expected outcome

A warning for non-monotonic data
A convenient option to restore the old behaviour, e.g. shading="drop"

cc @greglucas

Matplotlib version

Operating system: Linux Mint
Matplotlib version: '3.3.1'
Matplotlib backend (print(matplotlib.get_backend())): TkAgg
Python version: 3.8
Jupyter version (if applicable):
Other libraries:

I guess its possible matplotlib should be trying to do the grid interpolation in axes co-ordinates (and then back to data coordinates), not just data co-ordinates? I guess I'm a little ignorant as to how cartopy handles the wrap at 360 or +/-180. I assume they have a custom transform that handles the wrap? For matplotlib, it would be nice to make a scaled down version of a wrapping transform that we could include in our testing (since we won't test cartopy directly).

Of course, flat is just deprecated, but you may still use it. But under the auto kwarg, the convenient way to get the old behaviour is to do pcolormesh(x, y, Z[:-1, :-1]) where you can drop which column and row you would like. This was deemed preferable to magically dropping one of each with no warning or even documentation about which would be dropped. Better would be to do the interpolation yourself.

Of course, flat is just deprecated, but you may still use it.

For now. That's not a forward-looking work-around, though. The whole reason this comes up is because to avoid the deprecation warning, CartoPy is making 'auto' the default.

Cartopy generally doesn't handle the wrap. You need to plot in the right coordinate system or manually wrap your data. Really, it depends on the projection what the behavior really is around 0 or +/- 180.

I think the suggestion is that matplotlib might want to consider detecting when the 'nearest' adjustment is likely to produce artifacts like this, in cases where the coordinates are non-monotonic.

For now. That's not a forward-looking work-around, though.

Yes, the whole second paragraph was just to point out that the problem isn't urgent from Matplotlib's point of view and that there are workarounds.

FWIW, I'm not sure we have made anything worse or that anything has changed from Matplotlib's point of view. But I suggest we carry on the discussion in one spot, and since SciTools/cartopy#1638 has the most information, lets do it there?

I'll close this for now, just to keep the discussion in one spot, but not to claim there is no issue.

@dopplershift

I'll reopen as I think matplotlib should do something here. Not sure, however, if it is just warn/error on non-monotonic co-rodinates, which is what xarray does, or if we should try and interpolate in axes space rather than data space.

The example above and in SciTools/cartopy#1638 is a grid (or close to) in screen space, because the grid was defined on the same projection as the axes projection. So interpolating in screen space (i.e. transforming x and y finding the midpoints, and then transforming back) would be the right thing to do in that case. I think that would work in most cases.

@dopplershift, @QuLogic, @greglucas you all have more experience with non-linear projections than I do. Opinions? It seems xarray has struggled with this and come up with a set of levers via their infer_intervals kwarg. http://xarray.pydata.org/en/stable/generated/xarray.DataArray.plot.pcolormesh.html But in that case they know they may be on a map projection (which is the most reasonable use case for a non-linear projection?)

OK, here is a matplotlib-only repo using polar. This is like the error above, the grid is a cartesian grid, but specified in the polar co-ordinates. This is often done for geographic data sets where the units are specified in lat/lon, but the grid is really quadrilaterals in some projection.

import matplotlib.pyplot as plt
import numpy as np

x= np.arange(-1, 1.1, 0.4)
y= np.arange(-1, 1.1, 0.3)

X, Y = np.meshgrid(x, y)

pt = mproj.polar.PolarTransform()

R = X*0
T = X*0

for i in range(len(x)):
    for j in range(len(y)):
        R[j, i], T[j, i] = pt.inverted().transform((X[j, i], Y[j, i]))

ax = plt.subplot(121, projection='polar')
ax.plot(R, T, 'd')
ax.pcolormesh(R, T, R, shading='flat')

ax = plt.subplot(122, projection='polar')
ax.plot(R, T, 'd')
ax.pcolormesh(R, T, R, shading='nearest')

... another issue here is that our quadrilaterals are not centred very well for a non-regular grid. Probably should do something like a vectorized version of:

  newx[j, i] = (oldx[j, i], + oldx[j, i+1] + oldx[j+1, i] + oldx[j+1, i+1]) / 4
  newy[j, i] = (oldy[j, i], + oldy[j, i+1] + oldy[j+1, i] + oldy[j+1, i+1]) / 4

Here is a possible fix. It requires that the projection doesn't change int he life of the pcolormesh. I don't think there is a way to change an axes projection after it is created. (If anyone knows special numpy-fu to restagger the grid, feel free to comment. )

import matplotlib.pyplot as plt
import numpy as np

def make_grid(ax, X, Y):
    m, n = np.shape(X)
    
    # make the grid in display co-ordinates:
    Z = ax.transData.transform(np.vstack((X.flat, Y.flat)).T)
    X = np.reshape(Z[:, 0], (m,n))
    Y = np.reshape(Z[:, 1], (m, n))

    Xn = np.zeros((m+1, n+1))

    Yn = np.zeros((m+1, n+1))

    Xn[1:-1, 1:-1] = X[:-1, :-1]
    Xn[1:-1, 1:-1] += X[1:, :-1]
    Xn[1:-1, 1:-1] += X[1:, 1:]
    Xn[1:-1, 1:-1] += X[:-1, 1:]

    Xn[1:-1, 1:-1] *= 0.25
    # edges:
    Xn[:, 0] = 2 * Xn[:, 1] - Xn[:, 2]
    Xn[:, -1] = Xn[: , -2] + Xn[:, -2] - Xn[:, -3]   
    Xn[0, :] = Xn[1, :]
    Xn[-1, :] = Xn[-2, :]

    Yn[1:-1, 1:-1] = Y[:-1, :-1]
    Yn[1:-1, 1:-1] += Y[1:, :-1]
    Yn[1:-1, 1:-1] += Y[1:, 1:]
    Yn[1:-1, 1:-1] += Y[:-1, 1:]     
    Yn[1:-1, 1:-1] *= 0.25
    
    # edges:
    Yn[0, :] = 2* Yn[1, :] - Yn[2, :]
    Yn[-1, :] = 2* Yn[-2, :] - Yn[-3, :]
    Yn[:, 0] = Yn[:, 1]
    Yn[:, -1] = Yn[:, -2]
    
    # put back in original co-ordinates:
    Z = ax.transData.inverted().transform(np.vstack((Xn.flat, Yn.flat)).T)
    Xn = np.reshape(Z[:, 0], (m+1, n+1))
    Yn = np.reshape(Z[:, 1], (m+1, n+1))
   
    return Xn, Yn

    
fig, axs = plt.subplots(3, 1, figsize=(4,9), subplot_kw={'projection': 'polar'}, constrained_layout=True)
print(axs)
x= np.arange(0, 1.1, 0.4)
y= np.arange(-1, 1.1, 0.3)

X, Y = np.meshgrid(x, y)

pt = mproj.polar.PolarTransform()

ax=axs[0]
R = X * 0
T = X * 0

for i in range(len(x)):
    for j in range(len(y)):
        R[j, i], T[j, i] = pt.inverted().transform((X[j, i], Y[j, i]))

ax.plot(R, T, 'd')
ax.pcolormesh(R, T, R, shading='flat')
ax.set_title('Current "flat": drops a row')

ax=axs[1]

ax.plot(R, T, 'd')
ax.pcolormesh(R, T, R, shading='nearest')
ax.set_title('Current "nearest"')

Rn, Tn = make_grid(ax, R, T)
ax=axs[2]
ax.plot(R, T, 'd')
ax.pcolormesh(Rn, Tn, R, shading='flat')
ax.set_title('Proposed "nearest"')

I don't think this solution will work with Cartopy and the "wrapped" projection case? The original issue is that there is a 360 -> 0 "jump" in the middle of the x coordinates, so however you want to average these values, you're going to get something non-ideal because you actually want to account for the offset when calculating edges. If I have an array of coordinates that goes [357, 359, 1, 3] and those are actually the centers, then we'd need edges of [356, 358, 0, 2, 4] to make the mesh plot properly. This really boils down to originally defining pcolormesh with edges instead of centers.

I think we need to look into creating the mesh edge interpolation in Cartopy before sending anything to MPL. Right now, we are calling axes._pcolorargs('pcolormesh', ...) internally. I think the right solution may be to try and get away from using that private method and write our own logic instead.

I also think an MPL warning when coordinates are not strictly increasing/decreasing when you are creating new edges for people should probably be added.

It will work a wrapped projection if the data is monotonic in projection space (if the data is non-monotonic in projection space, flat will fail just as badly.)

If I have an array of coordinates that goes [357, 359, 1, 3] and those are actually the centers, then we'd need edges of [356, 358, 0, 2, 4]

Right, and thats what this proposed changed would give you. It would transform 357, 359, 1, 3 into screen units (say 100, 110, 120, 130) and find the mid points (say 95, 105, 115, 125, 135) and then transform those back into degrees (356, 358, 0, 2, 4).

Of course if the user gives data that transforms in screen space to 110, 120, 130, 100 then no method will work.

Ahh, you're right, I think that will fix the 360 -> 0 case because those are all increasing in screen space. But, for your last case, then in Cartopy we have the [179, 181, 183] case quite often. Those would map to [90, 10, 20] in screen space. We handle that by masking the wrapped coordinates, so we can't add an edge in there at 50.

We handle that by masking the wrapped coordinates, so we can't add an edge in there at 50.

and

Cartopy generally doesn't handle the wrap. You need to plot in the right coordinate system or manually wrap your data. Really, it depends on the projection what the behavior really is around 0 or +/- 180.

are contradictory. How do you handle the wrap 😉? It seems that if you masked the wrap before sending to pcolormesh, then the interpolation could put the masking in the right place to make things OK, but it'd be nice to have a simple example.

Here is a quick example to demonstrate both cases. Changing this to shading=flat looks acceptable in both cases, but drops the data.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import cartopy
import cartopy.crs as ccrs
print("MPL version: ", matplotlib.__version__)
print("Cartopy version: ", cartopy.__version__)

x1 = [179, 181, 183]
x2 = [359, 1, 3]
y = [-40, 0, 40]
xs1, ys = np.meshgrid(x1, y)
xs2, ys = np.meshgrid(x2, y)
z = np.arange(xs1.size).reshape(xs1.shape)

ax1 = plt.subplot(211, projection=ccrs.PlateCarree())
ax2 = plt.subplot(212, projection=ccrs.PlateCarree())

ax1.pcolormesh(xs1, ys, z, shading='auto', transform=ccrs.PlateCarree())
ax2.pcolormesh(xs2, ys, z, shading='auto', transform=ccrs.PlateCarree())

ax1.set_global()
ax2.set_global()

plt.show()

The wrap is handled by looking for "large" pcolormesh cells that extend over half the domain.
https://github.com/SciTools/cartopy/blob/8aca30fce1e9536aa002c8495995bccf375415a3/lib/cartopy/mpl/geoaxes.py#L1772-L1789
Those cells get masked and we then call pcolor with the original coordinates and Cartopy does some fancy boundary intersection to make one poly on each side of the boundary. Effectively turning the one quadmesh cell into two polys.

Phew, thats pretty complicated: ax2.plot(xs2, ys, 'o') won't plot the data point at 359 so far as I can tell. So its a bit hard to imagine the much more complicated pcolormesh working with a 359 in there...

So, it seems to me that the problem with my new algorithm and cartopy wraps is that the wrap masking logic is just not kicking in for small enough wraps. By interpolating the wraps are made smaller, i.e. 1/3 of the panel rather than 1/2, and they aren't being masked.

I don't think that is quite it. Even if we identified those cells as masked they would have a bad point in them and the pcolor algorithm would use that bad point to draw the polygons. I think that we need to address this on the Cartopy side because there are two competing data wraps (179 -> 181 and 359 -> 1). The way I think about it now is that Matplotlib is trying to be nice for users and make the edges that pcolormesh really expects. Cartopy in this case can be a user and calculate the edges before calling pcolormesh and negate the need to use the "auto" shading.

That would be ideal - or at worst case downstream libraries can drop the last row and column if they can't come up with a better way of specifying the corners.

OTOH, I think 1) matplotlib could still consider a screen space midpoint rather than a data space one. 2) matplotlib could warn if whatever interpolation we decide on is non-monotonic in any rows for X and any columns for Y

I believe this is closed by #18398....

thanks everyone!

Navigation Menu

Repository navigation

pcolormesh: shading='nearest' and non-monotonic coordinates #18317

Bug report

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

pcolormesh: shading='nearest' and non-monotonic coordinates #18317

Description

Bug report

Activity

jklymak commented on Aug 21, 2020

dopplershift commented on Aug 21, 2020

jklymak commented on Aug 21, 2020

jklymak commented on Aug 22, 2020

jklymak commented on Aug 22, 2020

jklymak commented on Aug 22, 2020

jklymak commented on Aug 23, 2020

greglucas commented on Aug 23, 2020

jklymak commented on Aug 23, 2020

greglucas commented on Aug 23, 2020

jklymak commented on Aug 23, 2020

greglucas commented on Aug 23, 2020

jklymak commented on Aug 23, 2020

jklymak commented on Aug 23, 2020

greglucas commented on Aug 23, 2020

jklymak commented on Aug 23, 2020

jklymak commented on Oct 28, 2020

mathause commented on Oct 28, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions