diff --git a/esmvaltool/diag_scripts/monitor/multi_datasets.py b/esmvaltool/diag_scripts/monitor/multi_datasets.py
index 068c4033da..70faee96c2 100644
--- a/esmvaltool/diag_scripts/monitor/multi_datasets.py
+++ b/esmvaltool/diag_scripts/monitor/multi_datasets.py
@@ -608,6 +608,7 @@
 from pprint import pformat
 
 import cartopy.crs as ccrs
+import dask.array as da
 import iris
 import matplotlib as mpl
 import matplotlib.dates as mdates
@@ -1178,8 +1179,15 @@ def _plot_map_with_ref(self, plot_func, dataset, ref_dataset):
             plot_kwargs['axes'] = axes_data
             if plot_func is iris.plot.contourf:
                 # see https://github.com/SciTools/cartopy/issues/2457
+                # and https://github.com/SciTools/cartopy/issues/2468
                 plot_kwargs['transform_first'] = True
-            plot_data = plot_func(cube, **plot_kwargs)
+                npx = da if cube.has_lazy_data() else np
+                cube_to_plot = cube.copy(
+                    npx.ma.filled(cube.core_data(), np.nan)
+                )
+            else:
+                cube_to_plot = cube
+            plot_data = plot_func(cube_to_plot, **plot_kwargs)
             axes_data.coastlines()
             if gridline_kwargs is not False:
                 axes_data.gridlines(**gridline_kwargs)
@@ -1196,7 +1204,17 @@ def _plot_map_with_ref(self, plot_func, dataset, ref_dataset):
             if self.plots[plot_type]['common_cbar']:
                 plot_kwargs.setdefault('vmin', plot_data.get_clim()[0])
                 plot_kwargs.setdefault('vmax', plot_data.get_clim()[1])
-            plot_ref = plot_func(ref_cube, **plot_kwargs)
+            if plot_func is iris.plot.contourf:
+                # see https://github.com/SciTools/cartopy/issues/2457
+                # and https://github.com/SciTools/cartopy/issues/2468
+                plot_kwargs['transform_first'] = True
+                npx = da if ref_cube.has_lazy_data() else np
+                ref_cube_to_plot = ref_cube.copy(
+                    npx.ma.filled(ref_cube.core_data(), np.nan)
+                )
+            else:
+                ref_cube_to_plot = ref_cube
+            plot_ref = plot_func(ref_cube_to_plot, **plot_kwargs)
             axes_ref.coastlines()
             if gridline_kwargs is not False:
                 axes_ref.gridlines(**gridline_kwargs)
@@ -1217,8 +1235,15 @@ def _plot_map_with_ref(self, plot_func, dataset, ref_dataset):
             plot_kwargs_bias['axes'] = axes_bias
             if plot_func is iris.plot.contourf:
                 # see https://github.com/SciTools/cartopy/issues/2457
+                # and https://github.com/SciTools/cartopy/issues/2468
                 plot_kwargs_bias['transform_first'] = True
-            plot_bias = plot_func(bias_cube, **plot_kwargs_bias)
+                npx = da if bias_cube.has_lazy_data() else np
+                bias_cube_to_plot = bias_cube.copy(
+                    npx.ma.filled(bias_cube.core_data(), np.nan)
+                )
+            else:
+                bias_cube_to_plot = bias_cube
+            plot_bias = plot_func(bias_cube_to_plot, **plot_kwargs_bias)
             axes_bias.coastlines()
             if gridline_kwargs is not False:
                 axes_bias.gridlines(**gridline_kwargs)
@@ -1276,8 +1301,15 @@ def _plot_map_without_ref(self, plot_func, dataset):
             plot_kwargs['axes'] = axes
             if plot_func is iris.plot.contourf:
                 # see https://github.com/SciTools/cartopy/issues/2457
+                # and https://github.com/SciTools/cartopy/issues/2468
                 plot_kwargs['transform_first'] = True
-            plot_map = plot_func(cube, **plot_kwargs)
+                npx = da if cube.has_lazy_data() else np
+                cube_to_plot = cube.copy(
+                    npx.ma.filled(cube.core_data(), np.nan)
+                )
+            else:
+                cube_to_plot = cube
+            plot_map = plot_func(cube_to_plot, **plot_kwargs)
             axes.coastlines()
             gridline_kwargs = self._get_gridline_kwargs(plot_type)
             if gridline_kwargs is not False:
diff --git a/esmvaltool/diag_scripts/shared/plot/_plot.py b/esmvaltool/diag_scripts/shared/plot/_plot.py
index 66f1e82c08..092479a999 100644
--- a/esmvaltool/diag_scripts/shared/plot/_plot.py
+++ b/esmvaltool/diag_scripts/shared/plot/_plot.py
@@ -4,6 +4,7 @@
 from copy import deepcopy
 
 import cartopy.crs as ccrs
+import dask.array as da
 import iris.quickplot
 import matplotlib.colors as colors
 import matplotlib.pyplot as plt
@@ -228,10 +229,17 @@ def global_contourf(cube,
     if cbar_range is not None:
         levels = np.linspace(*cbar_range)
         kwargs['levels'] = levels
-    kwargs['transform_first'] = True  # see SciTools/cartopy/issues/2457
     axes = plt.axes(projection=ccrs.Robinson(central_longitude=10))
     plt.sca(axes)
-    map_plot = iris.plot.contourf(cube, **kwargs)
+
+    # see https://github.com/SciTools/cartopy/issues/2457
+    # and https://github.com/SciTools/cartopy/issues/2468
+    kwargs['transform_first'] = True
+    npx = da if cube.has_lazy_data() else np
+    map_plot = iris.plot.contourf(
+        cube.copy(npx.ma.filled(cube.core_data(), np.nan)),
+        **kwargs,
+    )
 
     # Appearance
     axes.gridlines(color='lightgrey', alpha=0.5)