diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
index 17c553e6858..76cde15b9e1 100644
--- a/doc/changes/latest.inc
+++ b/doc/changes/latest.inc
@@ -109,6 +109,8 @@ Bugs
 
 - Fix bug with compensated CTF data when picking channels without preload (:gh:`8318` by `Eric Larson`_)
 
+- Fix bug with plotting MEG topographies where the wrong extrapolation made was used in ICA (:gh:`8637` by `Eric Larson`_)
+
 - Fix bug when merging fNIRS channels in :func:`mne.viz.plot_evoked_topomap` and related functions (:gh:`8306` by `Robert Luke`_)
 
 - Fix bug where events could overflow when writing to FIF (:gh:`8448` by `Eric Larson`_)
diff --git a/examples/preprocessing/plot_run_ica.py b/examples/preprocessing/plot_run_ica.py
index fe5efa07e03..ecd665fa670 100644
--- a/examples/preprocessing/plot_run_ica.py
+++ b/examples/preprocessing/plot_run_ica.py
@@ -10,7 +10,6 @@
 .. note:: This example does quite a bit of processing, so even on a
           fast machine it can take about a minute to complete.
 """
-
 # Authors: Denis Engemann <denis.engemann@gmail.com>
 #
 # License: BSD (3-clause)
diff --git a/examples/visualization/plot_evoked_topomap.py b/examples/visualization/plot_evoked_topomap.py
index 59185f92a9f..465e4ed282e 100644
--- a/examples/visualization/plot_evoked_topomap.py
+++ b/examples/visualization/plot_evoked_topomap.py
@@ -122,8 +122,8 @@
 # Animating the topomap
 # ---------------------
 #
-# Instead of using a still image we can plot magnetometer data as an animation
-# (animates only in matplotlib interactive mode)
+# Instead of using a still image we can plot magnetometer data as an animation,
+# which animates properly only in matplotlib interactive mode.
 
 # sphinx_gallery_thumbnail_number = 9
 times = np.arange(0.05, 0.151, 0.01)
diff --git a/mne/evoked.py b/mne/evoked.py
index c3527aba76e..f9181476cd6 100644
--- a/mne/evoked.py
+++ b/mne/evoked.py
@@ -381,7 +381,8 @@ def plot_joint(self, times="peaks", title='', picks=None,
     @fill_doc
     def animate_topomap(self, ch_type=None, times=None, frame_rate=None,
                         butterfly=False, blit=True, show=True, time_unit='s',
-                        sphere=None):
+                        sphere=None, *, extrapolate=_EXTRAPOLATE_DEFAULT,
+                        verbose=None):
         """Make animation of evoked data as topomap timeseries.
 
         The animation can be paused/resumed with left mouse button.
@@ -418,6 +419,10 @@ def animate_topomap(self, ch_type=None, times=None, frame_rate=None,
 
             .. versionadded:: 0.16
         %(topomap_sphere_auto)s
+        %(topomap_extrapolate)s
+
+            .. versionadded:: 0.22
+        %(verbose_meth)s
 
         Returns
         -------
@@ -433,7 +438,7 @@ def animate_topomap(self, ch_type=None, times=None, frame_rate=None,
         return _topomap_animation(
             self, ch_type=ch_type, times=times, frame_rate=frame_rate,
             butterfly=butterfly, blit=blit, show=show, time_unit=time_unit,
-            sphere=sphere)
+            sphere=sphere, extrapolate=extrapolate, verbose=verbose)
 
     def as_type(self, ch_type='grad', mode='fast'):
         """Compute virtual evoked using interpolated fields.
diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py
index df0ba4e935e..6871dbf08b9 100644
--- a/mne/preprocessing/ica.py
+++ b/mne/preprocessing/ica.py
@@ -1800,7 +1800,7 @@ def plot_components(self, picks=None, ch_type=None, res=64,
                         contours=6, image_interp='bilinear',
                         inst=None, plot_std=True, topomap_args=None,
                         image_args=None, psd_args=None, reject='auto',
-                        sphere=None):
+                        sphere=None, verbose=None):
         return plot_ica_components(self, picks=picks, ch_type=ch_type,
                                    res=res, vmin=vmin,
                                    vmax=vmax, cmap=cmap, sensors=sensors,
@@ -1810,19 +1810,21 @@ def plot_components(self, picks=None, ch_type=None, res=64,
                                    inst=inst, plot_std=plot_std,
                                    topomap_args=topomap_args,
                                    image_args=image_args, psd_args=psd_args,
-                                   reject=reject, sphere=sphere)
+                                   reject=reject, sphere=sphere,
+                                   verbose=verbose)
 
     @copy_function_doc_to_method_doc(plot_ica_properties)
     def plot_properties(self, inst, picks=None, axes=None, dB=True,
                         plot_std=True, topomap_args=None, image_args=None,
                         psd_args=None, figsize=None, show=True, reject='auto',
-                        reject_by_annotation=True):
+                        reject_by_annotation=True, *, verbose=None):
         return plot_ica_properties(self, inst, picks=picks, axes=axes,
                                    dB=dB, plot_std=plot_std,
                                    topomap_args=topomap_args,
                                    image_args=image_args, psd_args=psd_args,
                                    figsize=figsize, show=show, reject=reject,
-                                   reject_by_annotation=reject_by_annotation)
+                                   reject_by_annotation=reject_by_annotation,
+                                   verbose=verbose)
 
     @copy_function_doc_to_method_doc(plot_ica_sources)
     def plot_sources(self, inst, picks=None, start=None,
diff --git a/mne/viz/ica.py b/mne/viz/ica.py
index 852190bb745..d18ae389587 100644
--- a/mne/viz/ica.py
+++ b/mne/viz/ica.py
@@ -22,7 +22,7 @@
 from ..io.meas_info import create_info
 from ..io.pick import pick_types, _picks_to_idx
 from ..time_frequency.psd import psd_multitaper
-from ..utils import _reject_data_segments
+from ..utils import _reject_data_segments, verbose
 
 
 @fill_doc
@@ -250,11 +250,11 @@ def _get_psd_label_and_std(this_psd, dB, ica, num_std):
     return psd_ylabel, psds_mean, spectrum_std
 
 
-@fill_doc
+@verbose
 def plot_ica_properties(ica, inst, picks=None, axes=None, dB=True,
                         plot_std=True, topomap_args=None, image_args=None,
                         psd_args=None, figsize=None, show=True, reject='auto',
-                        reject_by_annotation=True):
+                        reject_by_annotation=True, *, verbose=None):
     """Display component properties.
 
     Properties include the topography, epochs image, ERP/ERF, power
@@ -308,6 +308,7 @@ def plot_ica_properties(ica, inst, picks=None, axes=None, dB=True,
     %(reject_by_annotation_raw)s
 
         .. versionadded:: 0.21.0
+    %(verbose)s
 
     Returns
     -------
diff --git a/mne/viz/tests/test_ica.py b/mne/viz/tests/test_ica.py
index b6831db68a0..0bc09886a8f 100644
--- a/mne/viz/tests/test_ica.py
+++ b/mne/viz/tests/test_ica.py
@@ -14,7 +14,7 @@
                  make_fixed_length_events)
 from mne.io import read_raw_fif
 from mne.preprocessing import ICA, create_ecg_epochs, create_eog_epochs
-from mne.utils import (run_tests_if_main, requires_sklearn, _click_ch_name,
+from mne.utils import (requires_sklearn, _click_ch_name, catch_logging,
                        _close_event)
 from mne.viz.ica import _create_properties_layout, plot_ica_properties
 from mne.viz.utils import _fake_click
@@ -70,7 +70,12 @@ def test_plot_ica_components():
     plt.close('all')
 
     # test interactive mode (passing 'inst' arg)
-    ica.plot_components([0, 1], image_interp='bilinear', inst=raw, res=16)
+    with catch_logging() as log:
+        ica.plot_components([0, 1], image_interp='bilinear', inst=raw, res=16,
+                            verbose='debug', ch_type='grad')
+    log = log.getvalue()
+    assert 'grad data' in log
+    assert 'Interpolation mode local to mean' in log
     fig = plt.gcf()
 
     # test title click
@@ -133,7 +138,11 @@ def test_plot_ica_properties():
         _create_properties_layout(figsize=(2, 2), fig=fig)
 
     topoargs = dict(topomap_args={'res': 4, 'contours': 0, "sensors": False})
-    ica.plot_properties(raw, picks=0, **topoargs)
+    with catch_logging() as log:
+        ica.plot_properties(raw, picks=0, verbose='debug', **topoargs)
+    log = log.getvalue()
+    assert raw.ch_names[0] == 'MEG 0113'
+    assert 'Interpolation mode local to mean' in log, log
     ica.plot_properties(epochs, picks=1, dB=False, plot_std=1.5, **topoargs)
     ica.plot_properties(epochs, picks=1, image_args={'sigma': 1.5},
                         topomap_args={'res': 4, 'colorbar': True},
@@ -396,6 +405,3 @@ def test_plot_instance_components():
     _fake_click(fig, ax, [line.get_xdata()[0], line.get_ydata()[0]], 'data')
     _fake_click(fig, ax, [-0.1, 0.9])  # click on y-label
     fig.canvas.key_press_event('escape')
-
-
-run_tests_if_main()
diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py
index 5b89ee64b79..6ee9b2af9fa 100644
--- a/mne/viz/tests/test_topomap.py
+++ b/mne/viz/tests/test_topomap.py
@@ -133,22 +133,27 @@ def test_plot_projs_topomap():
         plot_projs_topomap([eeg_proj], info_meg)
 
 
-def test_plot_topomap_animation():
+def test_plot_topomap_animation(capsys):
     """Test topomap plotting."""
     # evoked
     evoked = read_evokeds(evoked_fname, 'Left Auditory',
                           baseline=(None, 0))
     # Test animation
     _, anim = evoked.animate_topomap(ch_type='grad', times=[0, 0.1],
-                                     butterfly=False, time_unit='s')
+                                     butterfly=False, time_unit='s',
+                                     verbose='debug')
     anim._func(1)  # _animate has to be tested separately on 'Agg' backend.
+    out, _ = capsys.readouterr()
+    assert 'Interpolation mode local to 0' in out
     plt.close('all')
 
 
-def test_plot_topomap_animation_nirs(fnirs_evoked):
+def test_plot_topomap_animation_nirs(fnirs_evoked, capsys):
     """Test topomap plotting for nirs data."""
-    fig, anim = fnirs_evoked.animate_topomap(ch_type='hbo')
+    fig, anim = fnirs_evoked.animate_topomap(ch_type='hbo', verbose='debug')
     anim._func(1)  # _animate has to be tested separately on 'Agg' backend.
+    out, _ = capsys.readouterr()
+    assert 'Interpolation mode head to 0' in out
     assert len(fig.axes) == 2
     plt.close('all')
 
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index 9ecbf085608..828e4402ad4 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -767,6 +767,7 @@ def plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True,
 
 
 def _setup_interp(pos, res, extrapolate, sphere, outlines, border):
+    logger.debug(f'Interpolation mode {extrapolate} to {border}')
     xlim = np.inf, -np.inf,
     ylim = np.inf, -np.inf,
     mask_ = np.c_[outlines['mask_pos']]
@@ -792,6 +793,30 @@ def _setup_interp(pos, res, extrapolate, sphere, outlines, border):
     return extent, Xi, Yi, interp
 
 
+def _get_patch(outlines, extrapolate, interp, ax):
+    from matplotlib import patches
+    clip_radius = outlines['clip_radius']
+    clip_origin = outlines.get('clip_origin', (0., 0.))
+    _use_default_outlines = any(k.startswith('head') for k in outlines)
+    patch_ = None
+    if 'patch' in outlines:
+        patch_ = outlines['patch']
+        patch_ = patch_() if callable(patch_) else patch_
+        patch_.set_clip_on(False)
+        ax.add_patch(patch_)
+        ax.set_transform(ax.transAxes)
+        ax.set_clip_path(patch_)
+    if _use_default_outlines:
+        if extrapolate == 'local':
+            patch_ = patches.Polygon(
+                interp.mask_pts, clip_on=True, transform=ax.transData)
+        else:
+            patch_ = patches.Ellipse(
+                clip_origin, 2 * clip_radius[0], 2 * clip_radius[1],
+                clip_on=True, transform=ax.transData)
+    return patch_
+
+
 def _plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True,
                   res=64, axes=None, names=None, show_names=False, mask=None,
                   mask_params=None, outlines='head',
@@ -831,10 +856,7 @@ def _plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True,
             picks = list(range(data.shape[0]))
             pos = _find_topomap_coords(pos, picks=picks, sphere=sphere)
 
-    _check_option('extrapolate', extrapolate, ('box', 'local', 'head', 'auto'))
-    if extrapolate == 'auto':
-        extrapolate = 'local' if ch_type in _MEG_CH_TYPES_SPLIT else 'head'
-
+    extrapolate = _check_extrapolate(extrapolate, ch_type)
     if data.ndim > 1:
         raise ValueError("Data needs to be array of shape (n_sensors,); got "
                          "shape %s." % str(data.shape))
@@ -875,35 +897,16 @@ def _plot_topomap(data, pos, vmin=None, vmax=None, cmap=None, sensors=True,
     ax = axes if axes else plt.gca()
     _prepare_topomap(pos, ax)
 
-    _use_default_outlines = any(k.startswith('head') for k in outlines)
     mask_params = _handle_default('mask_params', mask_params)
 
     # find mask limits
-    clip_radius = outlines['clip_radius']
-    clip_origin = outlines.get('clip_origin', (0., 0.))
     extent, Xi, Yi, interp = _setup_interp(
         pos, res, extrapolate, sphere, outlines, border)
     interp.set_values(data)
     Zi = interp.set_locations(Xi, Yi)()
 
     # plot outline
-    patch_ = None
-    if 'patch' in outlines:
-        patch_ = outlines['patch']
-        patch_ = patch_() if callable(patch_) else patch_
-        patch_.set_clip_on(False)
-        ax.add_patch(patch_)
-        ax.set_transform(ax.transAxes)
-        ax.set_clip_path(patch_)
-    if _use_default_outlines:
-        from matplotlib import patches
-        if extrapolate == 'local':
-            patch_ = patches.Polygon(
-                interp.mask_pts, clip_on=True, transform=ax.transData)
-        else:
-            patch_ = patches.Ellipse(
-                clip_origin, 2 * clip_radius[0], 2 * clip_radius[1],
-                clip_on=True, transform=ax.transData)
+    patch_ = _get_patch(outlines, extrapolate, interp, ax)
 
     # plot interpolated map
     im = ax.imshow(Zi, cmap=cmap, vmin=vmin, vmax=vmax, origin='lower',
@@ -1010,7 +1013,7 @@ def _plot_ica_topomap(ica, idx=0, ch_type=None, res=64,
         data.ravel(), pos, vmin=vmin_, vmax=vmax_, res=res, axes=axes,
         cmap=cmap, outlines=outlines, contours=contours, sensors=sensors,
         image_interp=image_interp, show=show, extrapolate=extrapolate,
-        sphere=sphere, border=border)[0]
+        sphere=sphere, border=border, ch_type=ch_type)[0]
     if colorbar:
         cbar, cax = _add_colorbar(axes, im, cmap, pad=.05, title="AU",
                                   format='%3.2f')
@@ -1019,7 +1022,7 @@ def _plot_ica_topomap(ica, idx=0, ch_type=None, res=64,
     _hide_frame(axes)
 
 
-@fill_doc
+@verbose
 def plot_ica_components(ica, picks=None, ch_type=None, res=64,
                         vmin=None, vmax=None, cmap='RdBu_r',
                         sensors=True, colorbar=False, title=None,
@@ -1027,7 +1030,7 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64,
                         image_interp='bilinear',
                         inst=None, plot_std=True, topomap_args=None,
                         image_args=None, psd_args=None, reject='auto',
-                        sphere=None):
+                        sphere=None, *, verbose=None):
     """Project mixing matrix on interpolated sensor topography.
 
     Parameters
@@ -1109,6 +1112,7 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64,
         which applies the rejection parameters used when fitting
         the ICA object.
     %(topomap_sphere_auto)s
+    %(verbose)s
 
     Returns
     -------
@@ -1142,18 +1146,13 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64,
         figs = []
         for k in range(0, n_components, p):
             picks = range(k, min(k + p, n_components))
-            fig = plot_ica_components(ica, picks=picks, ch_type=ch_type,
-                                      res=res, vmax=vmax,
-                                      cmap=cmap, sensors=sensors,
-                                      colorbar=colorbar, title=title,
-                                      show=show, outlines=outlines,
-                                      contours=contours,
-                                      image_interp=image_interp, inst=inst,
-                                      plot_std=plot_std,
-                                      topomap_args=topomap_args,
-                                      image_args=image_args,
-                                      psd_args=psd_args, reject=reject,
-                                      sphere=sphere)
+            fig = plot_ica_components(
+                ica, picks=picks, ch_type=ch_type, res=res, vmax=vmax,
+                cmap=cmap, sensors=sensors, colorbar=colorbar, title=title,
+                show=show, outlines=outlines, contours=contours,
+                image_interp=image_interp, inst=inst, plot_std=plot_std,
+                topomap_args=topomap_args, image_args=image_args,
+                psd_args=psd_args, reject=reject, sphere=sphere)
             figs.append(fig)
         return figs
     else:
@@ -1164,7 +1163,7 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64,
     data = np.dot(ica.mixing_matrix_[:, picks].T,
                   ica.pca_components_[:ica.n_components_])
 
-    data_picks, pos, merge_channels, names, _, sphere, clip_origin = \
+    data_picks, pos, merge_channels, names, ch_type, sphere, clip_origin = \
         _prepare_topomap_plot(ica, ch_type, sphere=sphere)
     outlines = _make_head_outlines(sphere, pos, outlines, clip_origin)
 
@@ -1187,7 +1186,8 @@ def plot_ica_components(ica, picks=None, ch_type=None, res=64,
         im = plot_topomap(
             data_.flatten(), pos, vmin=vmin_, vmax=vmax_, res=res, axes=ax,
             cmap=cmap[0], outlines=outlines, contours=contours,
-            image_interp=image_interp, show=False, sensors=sensors)[0]
+            image_interp=image_interp, show=False, sensors=sensors,
+            ch_type=ch_type, **topomap_args)[0]
         im.axes.set_label(ica._ica_names[ii])
         if colorbar:
             cbar, cax = _add_colorbar(ax, im, cmap, title="AU",
@@ -2109,9 +2109,17 @@ def _hide_frame(ax):
     ax.set_frame_on(False)
 
 
-def _init_anim(ax, ax_line, ax_cbar, params, merge_channels, sphere):
+def _check_extrapolate(extrapolate, ch_type):
+    _check_option('extrapolate', extrapolate, ('box', 'local', 'head', 'auto'))
+    if extrapolate == 'auto':
+        extrapolate = 'local' if ch_type in _MEG_CH_TYPES_SPLIT else 'head'
+    return extrapolate
+
+
+@verbose
+def _init_anim(ax, ax_line, ax_cbar, params, merge_channels, sphere, ch_type,
+               extrapolate, verbose):
     """Initialize animated topomap."""
-    from matplotlib import pyplot as plt, patches
     logger.info('Initializing animation...')
     data = params['data']
     items = list()
@@ -2137,7 +2145,10 @@ def _init_anim(ax, ax_line, ax_cbar, params, merge_channels, sphere):
 
     _hide_frame(ax)
     extent, Xi, Yi, interp = _setup_interp(
-        params['pos'], 64, 'box', sphere, outlines, 0)
+        params['pos'], 64, extrapolate, sphere, outlines, 0)
+
+    patch_ = _get_patch(outlines, extrapolate, interp, ax)
+
     params['Zis'] = list()
     for frame in params['frames']:
         params['Zis'].append(interp.set_values(data[:, frame])(Xi, Yi))
@@ -2152,14 +2163,9 @@ def _init_anim(ax, ax_line, ax_cbar, params, merge_channels, sphere):
                    aspect='equal', extent=extent,
                    interpolation='bilinear')
     ax.autoscale(enable=True, tight=True)
-    plt.colorbar(im, cax=ax_cbar)
+    ax.figure.colorbar(im, cax=ax_cbar)
     cont = ax.contour(Xi, Yi, Zi, levels=cont_lims, colors='k', linewidths=1)
 
-    patch_ = patches.Ellipse((0, 0),
-                             2 * outlines['clip_radius'][0],
-                             2 * outlines['clip_radius'][1],
-                             clip_on=True,
-                             transform=ax.transData)
     im.set_clip_path(patch_)
     text = ax.text(0.55, 0.95, '', transform=ax.transAxes, va='center',
                    ha='right')
@@ -2249,7 +2255,7 @@ def _key_press(event, params):
 
 
 def _topomap_animation(evoked, ch_type, times, frame_rate, butterfly, blit,
-                       show, time_unit, sphere):
+                       show, time_unit, sphere, extrapolate, *, verbose=None):
     """Make animation of evoked data as topomap timeseries.
 
     See mne.evoked.Evoked.animate_topomap.
@@ -2273,7 +2279,6 @@ def _topomap_animation(evoked, ch_type, times, frame_rate, butterfly, blit,
         raise ValueError('All times must be inside the evoked time series.')
     frames = [np.abs(evoked.times - time).argmin() for time in times]
 
-    blit = False if plt.get_backend() == 'MacOSX' else blit
     picks, pos, merge_channels, _, ch_type, sphere, clip_origin = \
         _prepare_topomap_plot(evoked, ch_type, sphere=sphere)
     data = evoked.data[picks, :]
@@ -2292,6 +2297,7 @@ def _topomap_animation(evoked, ch_type, times, frame_rate, butterfly, blit,
         frames = np.linspace(0, len(evoked.times) - 1, frames).astype(int)
     ax_cbar = plt.subplot2grid(shape, (0, colspan), rowspan=rowspan)
     ax_cbar.set_title(_handle_default('units')[ch_type], fontsize=10)
+    extrapolate = _check_extrapolate(extrapolate, ch_type)
 
     params = dict(data=data, pos=pos, all_times=evoked.times, frame=0,
                   frames=frames, butterfly=butterfly, blit=blit,
@@ -2299,7 +2305,8 @@ def _topomap_animation(evoked, ch_type, times, frame_rate, butterfly, blit,
                   clip_origin=clip_origin)
     init_func = partial(_init_anim, ax=ax, ax_cbar=ax_cbar, ax_line=ax_line,
                         params=params, merge_channels=merge_channels,
-                        sphere=sphere)
+                        sphere=sphere, ch_type=ch_type,
+                        extrapolate=extrapolate, verbose=verbose)
     animate_func = partial(_animate, ax=ax, ax_line=ax_line, params=params)
     pause_func = partial(_pause_anim, params=params)
     fig.canvas.mpl_connect('button_press_event', pause_func)
diff --git a/tutorials/preprocessing/plot_40_artifact_correction_ica.py b/tutorials/preprocessing/plot_40_artifact_correction_ica.py
index c5b9f3c92d9..afa9af9f6dc 100644
--- a/tutorials/preprocessing/plot_40_artifact_correction_ica.py
+++ b/tutorials/preprocessing/plot_40_artifact_correction_ica.py
@@ -20,7 +20,6 @@
 and classes from that submodule:
 """
 
-
 import os
 import mne
 from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs,