pb12.py

ds_sigma_map = ds_sigma.sel(
  time='2000-01',
  sigma_0=26.6,
  method='nearest'
).sel(y=slice(-90,0)).squeeze()

import cartopy.crs as ccrs
import matplotlib.path as mpath
import cartopy.feature as cfeature

fig, axarr = plt.subplots(nrows=2,
                          figsize=[16,8],
                          subplot_kw={'projection':ccrs.PlateCarree(200)}
                         )
for ax, var, cmap in zip(
  axarr.flat,
  ['so', 'cfc11'],
  [cmo.haline, cmo.matter]
):
    
    if var == 'cfc11':
        vmin = 0
        vmax = 6e-9
    else:
        vmin = None
        vmax = None
    
    kwargs = dict(x='lon', y='lat', ax=ax, transform=ccrs.PlateCarree())
    
    ds_sigma_map[var].plot.contourf(
        levels=31, robust=True, cmap=cmap, vmin=vmin, vmax=vmax,  **kwargs
    )
ax.set_title('')
    ax.add_feature(cfeature.LAND, edgecolor='0.3')