Cartopy SouthPolarStereo projection makes (missing) data cover ocean and land features

Question

I have a plotting function that takes in input a series of parameters and plots a map. Data are in xarray data array format. The function is the following:

import xarray            as xr
import matplotlib.pyplot as plt
import cartopy.crs       as ccrs
import cartopy.feature   as cfeature
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib            import colors
import matplotlib

xticks_map = [-180, -120, -60, 0, 60, 120, 180]
yticks_map = [- 80, - 60, -40, -30]
# projection = ccrs.SouthPolarStereo(central_longitude = 0)
projection = ccrs.PlateCarree(central_longitude = 0)

def plotMap(da, levels = 11, vmin = None, vmax = None, name = "", title = "", SAVE = False, SHOW = True, 
            dpi = 300, figsize = (20,15), fontsize = 35, proj = ccrs.PlateCarree(), cmap = "viridis",
            shrink = .5, glLabels = {"top" : False, "bottom" : True, "right" : False, "left" : False}):
    matplotlib.rcParams.update({"font.size": fontsize})
    pt  = ccrs.PlateCarree()
    fig = plt.figure(dpi = dpi, figsize = figsize)
    ax  = plt.axes(  projection=proj)

    cbarKwg = {"shrink":shrink}
    da.plot.contourf(ax=ax, vmin = vmin, vmax = vmax, levels = levels, cmap = cmap, transform = ccrs.PlateCarree(), 
                          cbar_kwargs=cbarKwg)
    gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle = "--")
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.top_labels    = glLabels["top"]
    gl.bottom_labels = glLabels["bottom"]
    gl.right_labels  = glLabels["right"]
    gl.left_labels   = glLabels["left"]
    
    extent = {"xTicks" : xticks_map, "yTicks" : yticks_map}
    bounds = [np.min(extent["xTicks"]), np.max(extent["xTicks"]), np.min(extent["yTicks"]), np.max(extent["yTicks"])]
    ax.set_extent(bounds, crs = pt)
    gl.xlocator = mticker.FixedLocator(extent["xTicks"])
    gl.ylocator = mticker.FixedLocator(extent["yTicks"])
    
    ax.coastlines(resolution = "10m")
    ax.add_feature(cfeature.OCEAN)
    ax.add_feature(cfeature.LAND)
    ax.set_title(title)

    if SAVE:
        plt.tight_layout(pad=0)
        plt.savefig(name)
    if SHOW:
        plt.show()
    plt.close(fig)

My data are geographic data relative to a marine variable, hence over land there ara nan values. Same applies for latitudes outside the 80S - 30S belt. PlateCarree projection works just fine. When I use SouthPolarStereo (or Orthographic), where there should be nan, hence covered by the features OCEAN/LAND, there is a colour related to the 0 value in the colorbar. In those points there is actually the nan value (in fact PlateCarree works just fine), but somehow they are considered as normal values.

Now, the result with SouthPolarStereo projection is:

While with PlateCarree is:

Do you know a way out of this problem?

I noticed that if I coarsen the data: ds = ds.coarsen(dim={'latitude': 4, 'longitude': 4}, boundary='trim').mean() with a value of 4 and above, the result is the one I want. With 2 (and 1) I get the wrong output. Could it be that, since my data are high (not that much) resolution, the script get stuck somewhere? SouthPolarStereo projection takes a long time to process ( O(30 mins) ), while PlateCarree is much faster ( O(1 min) ). Data are like this:

>>> ds
<xarray.Dataset>
Dimensions:    (latitude: 200, longitude: 1799)
Coordinates:
  * latitude   (latitude) float64 -70.0 -69.8 -69.6 -69.4 ... -30.6 -30.4 -30.2
  * longitude  (longitude) float64 -179.8 -179.6 -179.4 ... 179.4 179.6 179.8
Data variables:
    FSLE       (latitude, longitude) float64 0.0 0.0 1.453 ... 0.0 7.763 nan

I am aware of this question But I could not use this way around in my case.