How to apply Earth Features and Land/Ocean masks with polar stereographic projections in Cartopy?

Question

I'm trying to mask my data to print a pic with polar stereographic projections

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')
#!MASK OCEAN
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()

but it didn't work ,it looks like the ocean cover the whole north hemisphere: enter image description here while if i try to mask the land withax.add_feature(cfeature.NaturalEarthFeature("physical", "land", "10m"), ec="red", fc="yellow", lw=2, alpha=0.4),it seems just right : enter image description here i'm so confused ,how can i mask the ocean properly in polar stereographic projections? is masking data the only way?

i tried to adjust the layer of ocean, finding the shape of ocean mask is not match. i used the same code but with different projection,

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()

it fuctioned. enter image description here

Answer 1

I find a solution after trying, first you need to get a shp file for ocean mask(the coastline file can't do):

import geopandas as gpd
sh1=gpd.read_file(r"F:\MASK\ne_10m_ocean\ne_10m_ocean.shp")
sh1_proj=sh1.to_crs(ccrs.NorthPolarStereo())
sh1_proj.plot(ax=ax, color='w',)

make this layer between your data and the gridline, that will make a satisfying mask. Maybe you don't want to set the zorder, because that may cause error on the display of gridline latitude label.

This solution doesn't fix the problem of cartopy, so it is still worth discussing.