pcolormesh shifts grid cells to the east incorrectly

[AgilentGCMS]

Description

I have a 2D array defined on a lat/lon grid that I can plot either with contourf or pcolormesh. If I use contourf the colored boxes show up exactly where expected. However, if I use pcolormesh the boxes are shifted to the east by a little bit. This is the same problem that's reported here. While that page offers a hacky fix, I'd much rather that cartopy did the right thing to begin with.

Code to reproduce

I am including the snippets I believe are relevant, but I can include a MWE if that helps. The 2D array to be plotted is called data, so first define a set of latitude and longitude array edges for pcolormesh

import numpy as np
nlat, nlon = data.shape
lat_edges = np.linspace(-90., 90., nlat+1)
lon_edges = np.linspace(-180., 180., nlon+1)

Now plot using pcolormesh (plot_ax has been set up before, as has data_proj and cmap)

vmin = np.percentile(data,  0.5)
vmax = np.percentile(data, 99.5)
img = plot_ax.pcolormesh(lon_edges, lat_edges, data, transform=data_proj, cmap=cmap, vmin=vmin, vmax=vmax, edgecolors='face', shading='flat')

This produces this map. There are two problems in this map:

1. The colored boxes are not bound by straight lat and lon lines, and as a result the map looks like an impressionist painting. This is more obvious if one zooms in.
2. The colored boxes are shifted a little bit to the east. This is obvious if one zooms in on the east and west coasts of (say) South America. In principle this could be a problem with the input data array, except that it's not, which will become clear when compared against contourf.

To plot with contourf, I do the following:

lon_centers = 0.5*(lon_edges[1:]+lon_edges[:-1])
lat_centers = 0.5*(lat_edges[1:]+lat_edges[:-1])
img = plot_ax.contourf(lon_centers, lat_centers, data, transform=data_proj, cmap=cmap, levels=np.linspace(vmin, vmax, 10), extend='both')

which produces this. Yes, there are some contouring artifacts, but neither of the above problems occur. Since my data are on a regular lat/lon grid, I'd really like to be able to just use pcolormesh and avoid contour calculations, but I can't because of the above problems.

Operating system

MacOS Ventura 13.4
(all python packages, including python itself, installed with macports)

Cartopy version

0.21.0

Matplotlib version

3.5.3

[greglucas]

Here is a MWE that works as I'd expect with no shifts, so I'm not following what you think needs updating.

import cartopy.crs as ccrs
from matplotlib import pyplot as plt
import numpy as np

# get the data
lats = [-10, 0, 10]
lons = [-10, 0, 10]

lon_center = 0.0
map_proj = ccrs.Robinson(central_longitude=lon_center)
data_proj = ccrs.PlateCarree(central_longitude=lon_center)

fig = plt.figure()
ax = plt.axes(projection=map_proj)
ax.coastlines()
# ax.set_global()
ax.gridlines(xlocs=lons, ylocs=lats, linewidth=0.5, linestyle=(0,(5,2)), color='0.25')
ax.pcolormesh(lons, lats, np.arange(4).reshape((2, 2)), transform=data_proj)

sc = ax.scatter(lons, lats, c='k', marker='o', transform=data_proj)

plt.show()

[rcomer]

This is the same problem that's reported here.

That page is saying that you need to specify lat/lon “edges” rather than “centres” to get what you want. You are already doing that. In any case that advice would now be out-of-date thanks to the shading options introduced in Matplotlib v3.3. So you should get something sensible regardless of whether you use “edges” or “centres”.
https://matplotlib.org/stable/users/prev_whats_new/whats_new_3.3.0.html#pcolor-and-pcolormesh-now-accept-shading-nearest-and-auto

[AgilentGCMS]

OK, I think I've figured it out. The problem reported here does not apply to be because I'm already specifying the lat/lon edges. Rather, my problem comes from how pcolormesh plots when I specify edgecolors='face'. Here's an MWE (ok, not really minimal).

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import pyplot as plt

# first, construct a global map
spherical_globe = ccrs.Globe(ellipse=None, semimajor_axis=6371000., semiminor_axis=6371000.)
map_proj = ccrs.Robinson(globe=spherical_globe, central_longitude=0.0)
data_proj = ccrs.PlateCarree(central_longitude=0.0)

aspect = (map_proj.x_limits[1]-map_proj.x_limits[0])/(map_proj.y_limits[1]-map_proj.y_limits[0])
height = 5.0
width = aspect * height
vpad = 0.3
cbar_pad = 0.5

figwidth = width
figheight = height + vpad + cbar_pad
fig = plt.figure(figsize=(figwidth, figheight))

plot_ax = plt.axes([0., (vpad+cbar_pad)/figheight, width/figwidth, height/figheight], projection=map_proj)
cbar_ax = plt.axes([0.05*width/figwidth, 0.6*cbar_pad/figheight, 0.9*width/figwidth, 0.4*cbar_pad/figheight])

plot_ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.7, edgecolor='k')
plot_ax.add_feature(cfeature.LAKES.with_scale('50m'), linewidth=0.5, edgecolor='k', facecolor='none')
plot_ax.gridlines(xlocs=np.arange(-180.,181.,30.), ylocs=np.arange(-90.,91.,30.), linewidth=0.5, linestyle=(0,(5,2)), color='0.25')

# now construct a "chessboard" data array
nlat = 360
nlon = 720
data = np.zeros((nlat,nlon), np.float32)
for i in range(nlat):
    for j in range(nlon):
        data[i,j] = (-1)**(i+j)

# now plot with pcolormesh
lat_edges = np.linspace(-90., 90., nlat+1)
lon_edges = np.linspace(-180., 180., nlon+1)

img = plot_ax.pcolormesh(lon_edges, lat_edges, data, transform=data_proj, cmap=plt.cm.bwr, vmin=-1.5, vmax=1.5, edgecolors='none', shading='flat')
cbar = plt.colorbar(mappable=img, cax=cbar_ax, orientation='horizontal')
plt.setp(cbar_ax.get_xticklabels(), size=16)

plt.show()

Now zoom into a small area at the intersection of (0,0), i.e., at the center of the plot, to minimize projection artifacts. Here is what edgecolors='none' looks like, while here is what edgecolors='face' looks like. Notice how, if I want edge colors to match the face colors, the blue squares bleed into blue squares to the NW and SE, while they're disconnected from the blue squares to the NE and SW? Same with the red squares. Furthermore, the boundary between the squares no longer aligns exactly with the lat/lon lines. This diagonal "bleeding" pattern is what made my pcolormesh plots look like impressionist paintings drawn with brush strokes going from the top left to the bottom right.

I don't know if there's a fix to this, or if a fix is even needed. But given that I've found a way around, which is to not draw the edges, I'm happy to close this issue.

[greglucas]

Just to add some more context: I don't think this is Cartopy specific, but can likely be traced up to Matplotlib. I'm guessing you're getting something similar to this issue: matplotlib/matplotlib#17323
and here is some longer commentary on the problem: matplotlib/matplotlib#1188