Search code examples
pythonmatplotlibcartopy

pcolormesh: artefacts / overlapping points when using ec='face'

I am plotting maps with pcolormesh and cartopy. When I don't use the option edgecolor or set it to None, results are like expected. But when I use ec='face', I get really strange artefacts and/or overlapping points - depending on the selected output format. I have tried with snap and antialiased but those don't help unfortunately. When I zoom into the interactive plot, it looks not that bad but my output file (pdf or png) looks like a mess.

How can I activate the edgecolor without having this weird behaviour?

MWE below. Using python 3.10.12, matplotlib 3.5.1, Cartopy 0.22.0.

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

proj1 = ccrs.PlateCarree()
proj2 = ccrs.LambertAzimuthalEqualArea()
fig   = plt.figure(figsize=(16,8), layout='constrained')

# generate coordinates
xx = np.linspace(-10,30,300)
yy = np.linspace( 35,75,300)
xxx,yyy = np.meshgrid(xx,yy)

rng   = np.random.default_rng(6666) 

# create artificial data
zzz = xxx + yyy 
noise = rng.normal(0,np.nanmax(zzz)*0.1,zzz.shape)
zzz += noise
 
# create some gaps in data
zdim = zzz.size
idx = rng.choice(range(zdim), size=zdim//3)

zshape  = zzz.shape
zz      = zzz.flatten()
zz[idx] = np.nan
zzz     = zz.reshape(zshape)

ax    = fig.add_subplot(1,2,1, projection=proj2)
ax.add_feature(cfeature.BORDERS,   lw=0.5, ec='k')
ax.add_feature(cfeature.COASTLINE, lw=0.5, ec='k')

# looks ok    
ax.pcolormesh(xxx, yyy, zzz,
              transform=proj1)

ax = fig.add_subplot(1,2,2, projection=proj2)
ax.add_feature(cfeature.BORDERS,   lw=0.5, ec='k')
ax.add_feature(cfeature.COASTLINE, lw=0.5, ec='k')

# looks not ok
ax.pcolormesh(xxx, yyy, zzz,
              #snap=True,
              #antialiased=True,
              ec='face',            # this seems to have a problem here
              transform=proj1)

plt.show()
#plt.savefig('plot_test.pdf')

enter image description here

Solution

  • The resolution of the data is quite high in comparison to the plot size/dpi. That combined with the fact that the edge is actually drawn at the edge, and not confined to the inside of a quadrilateral, causes the edges to overplot adjacent data.

    Increasing the figure size/dpi could help, or reducing the default linewidth of the edge to be more appropriate with respect to the size of the faces/quadrilaterals. Properties like linewidths are specified in points, which relate to the figure so it's not relative to the size/resolution of the data.

    Setting the linecolor to red emphasizes the behavior a bit better perhaps. The image below is with the grid reduces to 100x100, and zoomed in a little.

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

# generate coordinates
xx = np.linspace(-10,30, 100)
yy = np.linspace( 35,75, 100)
xxx,yyy = np.meshgrid(xx,yy)

rng   = np.random.default_rng(6666) 

# create artificial data
zzz = xxx + yyy 
noise = rng.normal(0,np.nanmax(zzz)*0.1,zzz.shape)
zzz += noise
 
# create some gaps in data
zdim = zzz.size
idx = rng.choice(range(zdim), size=zdim//3)

zshape  = zzz.shape
zz      = zzz.flatten()
zz[idx] = np.nan
zzz     = zz.reshape(zshape)

data_proj = ccrs.PlateCarree()
map_proj = ccrs.LambertAzimuthalEqualArea()

fig, axs = plt.subplot_mosaic(
    "ABC;DEF", figsize=(12,7), layout='constrained', 
    subplot_kw={"projection": map_proj}, facecolor="w",
)

axs["A"].set_title("No edge")
axs["A"].pcolormesh(xxx, yyy, zzz, transform=data_proj)

axs["B"].set_title("edge=face, default lw")
axs["B"].pcolormesh(xxx, yyy, zzz, transform=data_proj, ec='face')

axs["C"].set_title("edge=face, small lw")
axs["C"].pcolormesh(xxx, yyy, zzz, transform=data_proj, ec='face', lw=0.1)

axs["D"].set_title("No edge")
axs["D"].pcolormesh(xxx, yyy, zzz, transform=data_proj)

axs["E"].set_title("edge=r, default lw")
axs["E"].pcolormesh(xxx, yyy, zzz, transform=data_proj, ec='r')

axs["F"].set_title("edge=r, small lw")
axs["F"].pcolormesh(xxx, yyy, zzz, transform=data_proj, ec='r', lw=0.1)

for ax in axs.values():
    ax.add_feature(cfeature.BORDERS,   lw=0.5, ec='k', alpha=.5)
    ax.add_feature(cfeature.COASTLINE, lw=0.5, ec='k', alpha=.5)
    ax.set_extent((0, 30, 60, 76), crs=data_proj)

    enter image description here