set_boundary removes data from axis on faceting plots

I experienced that when using FacetGrid plots and if I want to apply a boundary to the axes, the data from the plots get removed.

A short example to illustrate this:

import xarray as xr
import matplotlib.path as mpath
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

ds = xr.tutorial.load_dataset('air_temperature')

data = ds.air.isel(time=range(0,2))
p = data.plot(transform=ccrs.PlateCarree(),
            y="lat", x="lon",
            col='time', col_wrap=2,
            subplot_kws={'projection': ccrs.NorthPolarStereo()}
)

theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)

for ax in p.axes.flat:
    ax.coastlines()
    ax.gridlines()
    ax.set_extent([-180,180, 50,90], crs=ccrs.PlateCarree())
    #ax.set_boundary(circle, transform=ax.transAxes)

With the last line I want to have a circle as a boundary of each axis. This works but the data are not shown anymore:

Commenting the last line in the example shows the data as expected but without the boundary of course:

Any hints on why this isn't working?

Cheers,
Markus

Start with the provided FacetGrid reproduction, especially the loop over p.axes.flat and the ax.set_boundary(circle, transform=ax.transAxes) call. Trace how the plotted data, Matplotlib axes, and Cartopy boundary interact; done means applying the circular boundary to every faceted axis without making the plotted data disappear.