在Cartopy北极极射投影中正确实现海洋掩膜的方法咨询

我最近在尝试用Cartopy绘制带北极极射投影的图像，想要对海洋区域做掩膜处理，但遇到了棘手的问题，想请教大家怎么解决。

先给大家看我最初用北极极射投影写的代码：

fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), subplot_kw={'projection': ccrs.NorthPolarStereo()}, dpi=300)
im = ax.pcolormesh(lon, lat, de, transform=ccrs.PlateCarree(), cmap='viridis')
# 添加国界
from cartopy.io.shapereader import Reader
shp = Reader('F:/MASK/global_all_country/global_all_country.shp')
ax.add_geometries(shp.geometries(), crs=ccrs.PlateCarree(), edgecolor='k', linewidths=0.5, facecolor='none')
ax.coastlines(resolution='10m')
# 尝试掩膜海洋
ax.add_feature(cfeature.NaturalEarthFeature("physical", "ocean", "50m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)
# 限制经纬度范围
box = [-180, 180, 15, 90]
ax.set_extent(box, crs=ccrs.PlateCarree())
# 添加圆形边界
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)
ax.set_boundary(circle, transform=ax.transAxes)
# 添加经纬度网格
gl = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), color='gray', linestyle='--', dms=True, linewidth=0.5,
                  xlocs=np.arange(-180, 180 + 60, 60), ylocs=np.arange(15, 90 + 15, 15))
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}
cbar = fig.colorbar(im, fraction=0.08, pad=0.1, location='bottom', aspect=30,
                        label=f'ΔAOD', ticks=[0.05, 0.1, 0.2, 0.3, 0.4, .5])
cbar.ax.tick_params(labelsize=10)
plt.show()

运行这段代码后，结果完全不对——海洋掩膜直接覆盖了整个北半球，根本没有正确对应实际的海洋区域。但奇怪的是，如果我把掩膜对象换成陆地，用下面这段代码：

ax.add_feature(cfeature.NaturalEarthFeature("physical", "land", "10m"), ec="red", fc="yellow", lw=2, alpha=0.4)

陆地的掩膜显示就完全正常，和实际陆地区域匹配。

我试过调整海洋图层的绘制顺序，但发现海洋掩膜的形状本身就和北极极射投影下的区域不匹配。不过当我把投影换成ccrs.PlateCarree()（正射投影）时，同样的海洋掩膜代码却能正常工作，代码如下：

fig, ax = plt.subplots(1, 1, figsize=(4, 3.5), subplot_kw={'projection': ccrs.PlateCarree()}, dpi=300)
im = ax.pcolormesh(lon, lat, de, transform=ccrs.PlateCarree(), cmap='viridis')
# 添加国界
from cartopy.io.shapereader import Reader
shp = Reader('F:/MASK/global_all_country/global_all_country.shp')
ax.add_geometries(shp.geometries(), crs=ccrs.PlateCarree(), edgecolor='k', linewidths=0.5, facecolor='none')
ax.coastlines(resolution='10m')
ax.add_feature(cfeature.NaturalEarthFeature("physical", "ocean", "50m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)
# 限制经纬度范围
box = [-180, 180, 15, 90]
ax.set_extent(box, crs=ccrs.PlateCarree())
# 添加经纬度网格
gl = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), color='gray', linestyle='--', dms=True, linewidth=0.5,
                  xlocs=np.arange(-180, 180 + 60, 60), ylocs=np.arange(15, 90 + 15, 15))
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10}
gl.ylabel_style = {'size': 10}
cbar = fig.colorbar(im, fraction=0.08, pad=0.1, location='bottom', aspect=30,
                        label=f'ΔAOD', ticks=[0.05, 0.1, 0.2, 0.3, 0.4, .5])
cbar.ax.tick_params(labelsize=10)
plt.show()

这个时候海洋掩膜就完全正常，只覆盖实际的海洋区域。

我现在特别困惑，想问问大家：在北极极射投影下，怎么才能正确实现海洋掩膜？难道只能通过直接对数据本身做掩膜的方法来实现吗？

备注：内容来源于stack exchange，提问作者rosco