diff --git a/doc/changes/dev/14220.newfeature.rst b/doc/changes/dev/14220.newfeature.rst
new file mode 100644
index 00000000000..73e00d98459
--- /dev/null
+++ b/doc/changes/dev/14220.newfeature.rst
@@ -0,0 +1 @@
+Add support for applying selected SSP projectors and reconstructing projected data in sensor space on :class:`~mne.io.Raw`, :class:`~mne.Epochs`, and :class:`~mne.Evoked` objects via ``apply_proj(projs=...)`` and the new public ``reconstruct_proj(projs=...)`` method, by `Hamza Abdelhedi`_.
diff --git a/mne/_fiff/proj.py b/mne/_fiff/proj.py
index e808f38de06..16418e2d601 100644
--- a/mne/_fiff/proj.py
+++ b/mne/_fiff/proj.py
@@ -3,6 +3,7 @@
 # Copyright the MNE-Python contributors.
 
 import re
+import warnings
 from copy import deepcopy
 from itertools import count
 
@@ -274,11 +275,14 @@ def add_proj(self, projs, remove_existing=False, verbose=None):
         return self
 
     @verbose
-    def apply_proj(self, verbose=None):
+    def apply_proj(self, *, projs=None, verbose=None):
         """Apply the signal space projection (SSP) operators to the data.
 
         Parameters
         ----------
+        projs : Projection | list of Projection | None
+            The projectors to apply. All projectors must already be present in
+            ``self.info["projs"]``. If ``None``, all projectors are applied.
         %(verbose)s
 
         Returns
@@ -288,9 +292,9 @@ def apply_proj(self, verbose=None):
 
         Notes
         -----
-        Once the projectors have been applied, they can no longer be
-        removed. It is usually not recommended to apply the projectors at
-        too early stages, as they are applied automatically later on
+        Once a projector has been applied, it can no longer be removed. It is
+        usually not recommended to apply the projectors at too early stages,
+        as they are applied automatically later on
         (e.g. when computing inverse solutions).
         Hint: using the copy method individual projection vectors
         can be tested without affecting the original data.
@@ -310,43 +314,85 @@ def apply_proj(self, verbose=None):
         from ..evoked import Evoked
         from ..io import BaseRaw
 
-        if self.info["projs"] is None or len(self.info["projs"]) == 0:
-            logger.info(
-                "No projector specified for this dataset. "
-                "Please consider the method self.add_proj."
-            )
-            return self
+        restore = None
+        if projs is not None:
+            if isinstance(projs, Projection):
+                projs = [projs]
+            projs = _check_projs(projs, copy=False)
+            if not projs:
+                return self
+            info = self.info.copy()
+            selected_idx = _proj_indices(info["projs"], projs)
+            to_apply = []
+            for ii in selected_idx:
+                proj = info["projs"][ii]
+                if proj["active"]:
+                    continue
+                # avoid emitting the warning twice
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore", RuntimeWarning)
+                    _, nproj, _ = make_projector([proj], info["ch_names"], info["bads"])
+                if nproj:
+                    to_apply.append(ii)
+            if not to_apply:
+                return self
+            with info._unlock():
+                attached = info["projs"]
+                keep = [
+                    proj["active"] or ii in to_apply for ii, proj in enumerate(attached)
+                ]
+                omitted = iter(proj for proj, use in zip(attached, keep) if not use)
+                info["projs"] = [proj for proj, use in zip(attached, keep) if use]
+            restore = (keep, omitted)
+            self.info = info
+
+        try:
+            if self.info["projs"] is None or len(self.info["projs"]) == 0:
+                logger.info(
+                    "No projector specified for this dataset. "
+                    "Please consider the method self.add_proj."
+                )
+                return self
 
-        # Exit delayed mode if you apply proj
-        if isinstance(self, BaseEpochs) and self._do_delayed_proj:
-            logger.info("Leaving delayed SSP mode.")
-            self._do_delayed_proj = False
+            # Exit delayed mode if you apply proj
+            if isinstance(self, BaseEpochs) and self._do_delayed_proj:
+                logger.info("Leaving delayed SSP mode.")
+                self._do_delayed_proj = False
 
-        if all(p["active"] for p in self.info["projs"]):
-            logger.info(
-                "Projections have already been applied. Setting proj attribute to True."
-            )
-            return self
+            if all(p["active"] for p in self.info["projs"]):
+                logger.info(
+                    "Projections have already been applied. "
+                    "Setting proj attribute to True."
+                )
+                return self
 
-        _projector, info = setup_proj(
-            deepcopy(self.info), add_eeg_ref=False, activate=True
-        )
-        # let's not raise a RuntimeError here, otherwise interactive plotting
-        if _projector is None:  # won't be fun.
-            logger.info("The projections don't apply to these data. Doing nothing.")
+            _projector, info = setup_proj(
+                self.info.copy(), add_eeg_ref=False, activate=True
+            )
+            # let's not raise a RuntimeError here, otherwise interactive plotting
+            if _projector is None:  # won't be fun.
+                logger.info("The projections don't apply to these data. Doing nothing.")
+                return self
+            self._projector, self.info = _projector, info
+            if isinstance(self, BaseRaw | Evoked):
+                if self.preload:
+                    self._data = np.dot(self._projector, self._data)
+            else:  # BaseEpochs
+                if self.preload:
+                    for ii, e in enumerate(self._data):
+                        self._data[ii] = self._project_epoch(e)
+                else:
+                    self.load_data()  # will automatically apply
+            logger.info("SSP projectors applied...")
             return self
-        self._projector, self.info = _projector, info
-        if isinstance(self, BaseRaw | Evoked):
-            if self.preload:
-                self._data = np.dot(self._projector, self._data)
-        else:  # BaseEpochs
-            if self.preload:
-                for ii, e in enumerate(self._data):
-                    self._data[ii] = self._project_epoch(e)
-            else:
-                self.load_data()  # will automatically apply
-        logger.info("SSP projectors applied...")
-        return self
+        finally:
+            if restore is not None:
+                keep, omitted = restore
+                visible = iter(self.info["projs"])
+                with self.info._unlock():
+                    self.info["projs"] = [
+                        next(visible) if use else next(omitted) for use in keep
+                    ]
 
     def del_proj(self, idx="all"):
         """Remove SSP projection vector.
@@ -515,19 +561,63 @@ def plot_projs_topomap(
         )
         return fig
 
-    def _reconstruct_proj(self, mode="accurate", origin="auto"):
+    def reconstruct_proj(self, *, projs=None, mode="accurate", origin="auto"):
+        """Apply SSP projectors and reconstruct the resulting signal in sensor space.
+
+        Operates in place.
+
+        Parameters
+        ----------
+        projs : Projection | list of Projection | None
+            The projector or projectors to apply before reconstruction. All
+            projectors must already be present in ``self.info["projs"]``. If
+            ``None``, all projectors attached to the instance are used.
+        mode : str
+            Either ``'accurate'`` or ``'fast'``, determines the quality of the
+            Legendre polynomial expansion used for reconstruction.
+        origin : array-like, shape (3,) | str
+            Origin of the sphere in the head coordinate frame and in meters.
+            Can be ``'auto'`` (default), which means a head-digitization-based
+            origin fit.
+
+        Returns
+        -------
+        self : same type as the input data
+            The modified instance.
+        """
         from ..forward import _map_meg_or_eeg_channels
 
-        if len(self.info["projs"]) == 0:
-            return self
-        self.apply_proj()
+        if projs is None:
+            if len(self.info["projs"]) == 0:
+                return self
+            self.apply_proj()
+            mapping_info = self.info
+            selected_projs = None
+        else:
+            self.apply_proj(projs=projs)
+            selected_projs = [projs] if isinstance(projs, Projection) else projs
+            if len(selected_projs) == 0:
+                return self
+            mapping_info = self.info.copy()
+            with mapping_info._unlock():
+                mapping_info["projs"] = [
+                    proj for proj in mapping_info["projs"] if proj["active"]
+                ]
         for kind in ("meg", "eeg"):
             kwargs = dict(meg=False)
             kwargs[kind] = True
             picks = pick_types(self.info, **kwargs)
             if len(picks) == 0:
                 continue
-            info_from = pick_info(self.info, picks)
+            info_from = pick_info(mapping_info, picks)
+            if selected_projs is not None:
+                with warnings.catch_warnings():
+                    warnings.simplefilter("ignore", RuntimeWarning)
+                    _, nproj, _ = make_projector(
+                        selected_projs, info_from["ch_names"], info_from["bads"]
+                    )
+                if nproj == 0:
+                    continue
             info_to = info_from.copy()
             with info_to._unlock():
                 info_to["projs"] = []
@@ -1164,6 +1254,27 @@ def setup_proj(
     return projector, info
 
 
+def _proj_indices(attached, projs):
+    """Find the indices of projectors in an attached projector list."""
+    selected = []
+    for pi, proj in enumerate(projs):
+        matches = [
+            ii
+            for ii, attached_proj in enumerate(attached)
+            if _proj_equal(proj, attached_proj, check_active=False)
+        ]
+        if len(matches) == 0:
+            raise ValueError(
+                f"projs[{pi}] does not match any projector in self.info['projs']"
+            )
+        if len(matches) > 1:
+            raise ValueError(
+                f"projs[{pi}] matches multiple projectors in self.info['projs']"
+            )
+        selected.append(matches[0])
+    return list(dict.fromkeys(selected))
+
+
 def _uniquify_projs(projs, check_active=True, sort=True):
     """Make unique projs."""
     final_projs = []
diff --git a/mne/tests/test_proj.py b/mne/tests/test_proj.py
index 8ed5b74c87a..316258c6eb0 100644
--- a/mne/tests/test_proj.py
+++ b/mne/tests/test_proj.py
@@ -12,6 +12,7 @@
 
 from mne import (
     Epochs,
+    EvokedArray,
     compute_proj_epochs,
     compute_proj_evoked,
     compute_proj_raw,
@@ -26,6 +27,7 @@
 )
 from mne._fiff.proj import (
     _EEG_AVREF_PICK_DICT,
+    Projection,
     _needs_eeg_average_ref_proj,
     activate_proj,
     make_projector,
@@ -57,6 +59,265 @@
 ecg_fname = sample_path / "sample_audvis_ecg-proj.fif"
 
 
+def _make_test_proj(ch_names, vector, desc):
+    """Make a projector for selection tests."""
+    vector = np.atleast_2d(np.asarray(vector, float))
+    return Projection(
+        data=dict(
+            col_names=ch_names,
+            row_names=None,
+            nrow=len(vector),
+            ncol=len(ch_names),
+            data=vector,
+        ),
+        desc=desc,
+    )
+
+
+def _make_selection_raw():
+    """Make mixed EEG/MEG data with three attached projectors."""
+    ch_names = [f"EEG {ii:03d}" for ii in range(3)]
+    ch_names += [f"MEG {ii:03d}" for ii in range(3)]
+    info = create_info(ch_names, 100.0, ["eeg"] * 3 + ["mag"] * 3)
+    data = np.random.default_rng(0).standard_normal((6, 600))
+    raw = RawArray(data, info, verbose=False)
+    projs = [
+        _make_test_proj(ch_names[:3], [1.0, 1.0, 0.0], "EEG-a"),
+        _make_test_proj(ch_names[:3], [0.0, 1.0, 1.0], "EEG-b"),
+        _make_test_proj(ch_names[3:], [1.0, 1.0, 0.0], "MEG"),
+    ]
+    raw.add_proj(projs, verbose=False)
+    return raw, projs
+
+
+def _active_projs(inst):
+    """Return the active state of attached projectors."""
+    return [proj["active"] for proj in inst.info["projs"]]
+
+
+def _make_selection_epochs(raw, *, preload=False, proj=False, reject=None, events=None):
+    """Make Epochs for projector selection tests."""
+    if events is None:
+        events = (100, 300, 500)
+    events = np.column_stack(
+        [events, np.zeros(len(events), int), np.ones(len(events), int)]
+    )
+    return Epochs(
+        raw,
+        events,
+        1,
+        0.0,
+        0.2,
+        baseline=None,
+        reject=reject,
+        proj=proj,
+        preload=preload,
+        verbose=False,
+    )
+
+
+def _make_selection_instance(kind):
+    """Make a Raw, Epochs, or Evoked instance for projection selection tests."""
+    raw, projs = _make_selection_raw()
+    if kind == "epochs":
+        inst = _make_selection_epochs(raw, preload=True)
+    elif kind == "evoked":
+        inst = EvokedArray(raw.get_data()[:, :10], raw.info, tmin=0.0)
+    else:
+        inst = raw
+    return inst, projs
+
+
+def _get_reconstruction_evoked(raw_orig, events, picks):
+    """Get an Evoked instance for reconstruction tests."""
+    raw = raw_orig.copy()
+    raw.add_proj([], remove_existing=True)
+    epochs = Epochs(
+        raw,
+        events[:5],
+        1,
+        -0.1,
+        0.1,
+        picks=picks,
+        decim=10,
+        verbose="error",
+    )
+    epochs.info["bads"] = [epochs.ch_names[-5], epochs.ch_names[-1]]
+    epochs.info.normalize_proj()
+    return epochs.average()
+
+
+def test_reconstruct_proj(raw_orig, events):
+    """Test default and selected SSP projector reconstruction."""
+    cases = (
+        (
+            (0, 1, 2, 3, 4, 6, 7, 61, 122, 183, 244, 305),
+            (0.63, 0.65),
+            False,
+        ),
+        (np.arange(340, 360), (0.56, 0.57), True),
+        (np.arange(340, 360), (0.79, 0.81), False),
+    )
+    for picks, rlims, avg_proj in cases:
+        evoked = _get_reconstruction_evoked(raw_orig, events, picks)
+        if avg_proj:
+            evoked.set_eeg_reference(projection=True).apply_proj(verbose=False)
+        original = evoked.data.copy()
+        if not avg_proj:
+            assert_allclose(evoked.copy().reconstruct_proj().data, original)
+
+        proj = compute_proj_evoked(evoked.copy().crop(None, 0).apply_proj())
+        evoked.add_proj(proj, verbose=False)
+        projected = evoked.copy().apply_proj(verbose=False).data
+        reconstructed = evoked.copy().reconstruct_proj().data
+        norm = np.linalg.norm(original)
+        norm_proj = np.linalg.norm(projected)
+        norm_recon = np.linalg.norm(reconstructed)
+        r = np.dot(reconstructed.ravel(), original.ravel()) / (norm_recon * norm)
+        assert rlims[0] < r < rlims[1]
+        assert 1.05 * norm_proj < norm_recon
+        if not avg_proj:
+            assert norm_proj < norm * 0.9
+
+    evoked = _get_reconstruction_evoked(
+        raw_orig,
+        events,
+        (0, 1, 2, 3, 4, 6, 7, 61, 122, 183, 244, 305, 315, 316, 317, 318),
+    )
+    original = evoked.data.copy()
+    meg_proj = compute_proj_evoked(
+        evoked.copy().pick("meg"), n_grad=0, n_mag=1, n_eeg=0
+    )[0]
+    eeg_proj = compute_proj_evoked(
+        evoked.copy().pick("eeg"), n_grad=0, n_mag=0, n_eeg=1
+    )[0]
+    evoked.add_proj([meg_proj, eeg_proj], verbose=False)
+    meg_picks = pick_types(evoked.info, meg=True, eeg=False)
+    eeg_picks = pick_types(evoked.info, meg=False, eeg=True)
+    reconstructed_all = evoked.copy().reconstruct_proj().data
+    for selected, selected_picks, untouched_picks in (
+        (meg_proj, meg_picks, eeg_picks),
+        (eeg_proj, eeg_picks, meg_picks),
+    ):
+        reconstructed = evoked.copy().reconstruct_proj(projs=selected).data
+        assert_allclose(reconstructed[untouched_picks], original[untouched_picks])
+        assert_allclose(
+            reconstructed[selected_picks], reconstructed_all[selected_picks]
+        )
+        assert not np.array_equal(
+            reconstructed_all[untouched_picks], original[untouched_picks]
+        )
+
+
+@pytest.mark.parametrize("kind", ["raw", "epochs", "evoked"])
+def test_apply_proj_default(kind):
+    """Test that ``projs=None`` preserves legacy behavior."""
+    inst, _ = _make_selection_instance(kind)
+    legacy = inst.copy().apply_proj(verbose=False)
+    explicit_none = inst.copy().apply_proj(projs=None, verbose=False)
+    assert_allclose(explicit_none.get_data(), legacy.get_data())
+    assert _active_projs(explicit_none) == [True] * 3
+
+
+@pytest.mark.parametrize("kind", ["raw", "epochs", "evoked"])
+@pytest.mark.parametrize("selected_idx", [(0,), (1,), (0, 1), (2,)])
+def test_apply_proj_selection(kind, selected_idx):
+    """Test one, multiple, and modality-specific projector selections."""
+    inst, projs = _make_selection_instance(kind)
+    data = inst.get_data().copy()
+    selected = [projs[ii] for ii in selected_idx]
+    passed = cp.deepcopy(selected[0] if len(selected) == 1 else selected)
+    got = inst.copy().apply_proj(projs=passed, verbose=False)
+    projector = make_projector(selected, inst.ch_names)[0]
+    assert_allclose(got.get_data(), np.matmul(projector, data))
+    assert_allclose(got._projector, projector)
+    assert _active_projs(got) == [ii in selected_idx for ii in range(3)]
+    assert not any(proj["active"] for proj in np.atleast_1d(passed))
+    if 2 in selected_idx:
+        assert_allclose(got.get_data()[..., :3, :], data[..., :3, :])
+    else:
+        assert_allclose(got.get_data()[..., 3:, :], data[..., 3:, :])
+    if kind == "evoked" and selected_idx == (0,):
+        expected = inst.copy().del_proj()
+        expected.add_proj(projs[0], verbose=False).apply_proj(verbose=False)
+        expected.add_proj(projs[1:], verbose=False)
+        assert_allclose(got._data, expected._data)
+        assert_allclose(got._projector, expected._projector)
+        assert got.info["projs"] == expected.info["projs"]
+
+
+def test_apply_proj_selection_invalid():
+    """Test non-attached and ambiguous projector selections."""
+    raw, projs = _make_selection_raw()
+    unattached = _make_test_proj(raw.ch_names[:3], [1.0, 0.0, 1.0], "other")
+    with pytest.raises(ValueError, match="does not match"):
+        raw.apply_proj(projs=unattached, verbose=False)
+
+    raw.info["projs"].append(cp.deepcopy(projs[0]))
+    with pytest.raises(ValueError, match="matches multiple"):
+        raw.apply_proj(projs=projs[0], verbose=False)
+
+
+def test_apply_proj_selection_restores_on_error():
+    """Test that omitted projectors are restored if projection fails."""
+    ch_names = ["EEG 001", "EEG 002"]
+    raw = RawArray(np.zeros((2, 10)), create_info(ch_names, 100.0, "eeg"))
+    projs = [
+        _make_test_proj(ch_names, [1.0, 0.0], "A"),
+        _make_test_proj(ch_names, [0.0, 1.0], "B"),
+        _make_test_proj(ch_names, [1.0, 1.0], "C"),
+    ]
+    raw.add_proj(projs, verbose=False).apply_proj(projs=projs[0], verbose=False)
+    with pytest.raises(RuntimeError, match="will yield no components"):
+        raw.apply_proj(projs=projs[1], verbose=False)
+    assert _active_projs(raw) == [True, False, False]
+    assert [proj["desc"] for proj in raw.info["projs"]] == ["A", "B", "C"]
+
+
+@pytest.mark.parametrize("kind", ["raw", "epochs", "evoked"])
+@pytest.mark.parametrize("selection", ["unsupported", "applicable"])
+def test_apply_proj_selection_after_active(kind, selection):
+    """Test selected projectors after an applicable projector is active."""
+    inst, projs = _make_selection_instance(kind)
+    original = inst.get_data().copy()
+    inst.apply_proj(projs=projs[0], verbose=False)
+    after_active = inst.get_data().copy()
+    active_projector = inst._projector.copy()
+
+    if selection == "unsupported":
+        selected = _make_test_proj(["missing"], [1.0], "unsupported")
+        inst.add_proj(selected, verbose=False)
+    else:
+        selected = projs[1]
+    inst.apply_proj(projs=selected, verbose=False)
+
+    if selection == "unsupported":
+        assert_allclose(inst.get_data(), after_active)
+        assert_allclose(inst._projector, active_projector)
+        assert _active_projs(inst) == [True, False, False, False]
+    else:
+        projector = make_projector(projs[:2], inst.ch_names)[0]
+        assert_allclose(inst.get_data(), np.matmul(projector, original))
+        assert_allclose(inst._projector, projector)
+        assert _active_projs(inst) == [True, True, False]
+
+
+@pytest.mark.parametrize("preload", [False, True])
+@pytest.mark.parametrize("proj", [False, "delayed"])
+def test_apply_proj_selection_epochs(preload, proj):
+    """Test selected projection for lazy, preloaded, and delayed Epochs."""
+    raw, projs = _make_selection_raw()
+    epochs = _make_selection_epochs(raw, preload=preload, proj=proj)
+    data = epochs.get_data(copy=True)
+    matrix = make_projector([projs[0]], epochs.ch_names)[0]
+    epochs.apply_proj(projs=projs[0], verbose=False)
+    assert epochs.preload
+    assert not epochs._do_delayed_proj
+    assert_allclose(epochs.get_data(copy=True), np.matmul(matrix, data))
+    assert_allclose(epochs._projector, matrix)
+    assert _active_projs(epochs) == [True, False, False]
+
+
 def test_bad_proj():
     """Test dealing with bad projection application."""
     raw = read_raw_fif(raw_fname, preload=True)
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 324d7f79baa..cee4e4d315e 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -450,7 +450,7 @@ def _plot_evoked(
     if projector is not None:
         evoked.data[:] = np.dot(projector, evoked.data)
     if proj == "reconstruct":
-        evoked = evoked._reconstruct_proj()
+        evoked = evoked.reconstruct_proj()
 
     if plot_type == "butterfly":
         _plot_lines(
@@ -1961,7 +1961,7 @@ def plot_evoked_joint(
     if proj:
         evoked.apply_proj()
         if proj == "reconstruct":
-            evoked._reconstruct_proj()
+            evoked.reconstruct_proj()
     topomap_args["proj"] = ts_args["proj"] = False  # don't reapply
     evoked.pick(picks, exclude=exclude)
     info = evoked.info
diff --git a/mne/viz/tests/test_evoked.py b/mne/viz/tests/test_evoked.py
index d1c37dacdb9..cd489d88416 100644
--- a/mne/viz/tests/test_evoked.py
+++ b/mne/viz/tests/test_evoked.py
@@ -17,7 +17,6 @@
 from mne import (
     Epochs,
     compute_covariance,
-    compute_proj_evoked,
     compute_rank,
     make_fixed_length_events,
     read_cov,
@@ -301,51 +300,14 @@ def _get_amplitudes(fig):
     return amplitudes
 
 
-@pytest.mark.parametrize(
-    "picks, rlims, avg_proj",
-    [
-        (default_picks[:-4], (0.59, 0.61), False),  # MEG
-        (np.arange(340, 360), (0.56, 0.57), True),  # EEG
-        (np.arange(340, 360), (0.79, 0.81), False),  # EEG
-    ],
-)
-def test_plot_evoked_reconstruct(picks, rlims, avg_proj):
+def test_plot_evoked_reconstruct():
     """Test proj="reconstruct"."""
-    evoked = _get_epochs(picks=picks).average()
-    if avg_proj:
-        evoked.set_eeg_reference(projection=True).apply_proj()
-        assert len(evoked.info["projs"]) == 1
-        assert evoked.proj is True
-    else:
-        assert len(evoked.info["projs"]) == 0
-        assert evoked.proj is False
-    fig = evoked.plot(
-        proj=True, hline=[1], exclude=[], window_title="foo", time_unit="s"
-    )
-    amplitudes = _get_amplitudes(fig)
-    assert len(amplitudes) == len(picks)
-    assert evoked.proj is avg_proj
+    evoked = _get_epochs(picks=np.arange(340, 360)).average()
+    evoked.set_eeg_reference(projection=True).apply_proj()
 
     fig = evoked.plot(proj="reconstruct", exclude=[])
     amplitudes_recon = _get_amplitudes(fig)
-    if avg_proj is False:
-        assert_allclose(amplitudes, amplitudes_recon)
-    proj = compute_proj_evoked(evoked.copy().crop(None, 0).apply_proj())
-    evoked.add_proj(proj)
-    assert len(evoked.info["projs"]) == 2 if len(picks) == 3 else 4
-    fig = evoked.plot(proj=True, exclude=[])
-    amplitudes_proj = _get_amplitudes(fig)
-    fig = evoked.plot(proj="reconstruct", exclude=[])
-    amplitudes_recon = _get_amplitudes(fig)
-    assert len(amplitudes_recon) == len(picks)
-    norm = np.linalg.norm(amplitudes)
-    norm_proj = np.linalg.norm(amplitudes_proj)
-    norm_recon = np.linalg.norm(amplitudes_recon)
-    r = np.dot(amplitudes_recon.ravel(), amplitudes.ravel()) / (norm_recon * norm)
-    assert rlims[0] < r < rlims[1]
-    assert 1.05 * norm_proj < norm_recon
-    if not avg_proj:
-        assert norm_proj < norm * 0.9
+    assert len(amplitudes_recon) == len(evoked.ch_names)
 
     cov = read_cov(cov_fname)
     with pytest.raises(ValueError, match='Cannot use proj="reconstruct"'):
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index 75ca4371cb6..3cad9b9b9b5 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -2528,7 +2528,7 @@ def _plot_evoked_topomap(
     if proj is True and not evoked.proj:
         evoked.apply_proj()
     elif proj == "reconstruct":
-        evoked._reconstruct_proj()
+        evoked.reconstruct_proj()
 
     # remove compensation matrices (safe: only plotting & already made copy)
     with evoked.info._unlock():
diff --git a/tutorials/preprocessing/45_projectors_background.py b/tutorials/preprocessing/45_projectors_background.py
index 3c83d49d8c3..392c44dc675 100644
--- a/tutorials/preprocessing/45_projectors_background.py
+++ b/tutorials/preprocessing/45_projectors_background.py
@@ -467,11 +467,16 @@ def setup_3d_axes():
 # automatically applied. It is also possible to apply projectors manually when
 # working with :class:`~mne.io.Raw`, :class:`~mne.Epochs` or
 # :class:`~mne.Evoked` objects via the object's :meth:`~mne.io.Raw.apply_proj`
-# method. For all instance types, you can always copy the contents of
-# :samp:`{<instance>}.info['projs']` into a separate :class:`list` variable,
-# use :samp:`{<instance>}.del_proj({<index of proj(s) to remove>})` to remove
-# one or more projectors, and then add them back later with
-# :samp:`{<instance>}.add_proj({<list containing projs>})` if desired.
+# method. By default, ``apply_proj()`` applies all projectors attached to the
+# instance. To apply only a particular projector or subset, pass the desired
+# projectors using ``projs=``. The selected projectors must already be attached
+# to the instance. For example, ``ecg_projs[2:4]`` are the two magnetometer
+# ("axial") ECG projectors, so we can apply only those:
+
+raw_ecg_mag = raw.copy().apply_proj(projs=ecg_projs[2:4])
+
+# %%
+# This avoids having to temporarily remove the other attached projectors.
 #
 # .. warning::
 #
diff --git a/tutorials/preprocessing/50_artifact_correction_ssp.py b/tutorials/preprocessing/50_artifact_correction_ssp.py
index e99a1dab7ff..8d1266c60a3 100644
--- a/tutorials/preprocessing/50_artifact_correction_ssp.py
+++ b/tutorials/preprocessing/50_artifact_correction_ssp.py
@@ -534,6 +534,18 @@
     for text in list(ax.texts):
         text.remove()
 
+# %%
+# Using ``proj="reconstruct"`` is convenient when reconstruction is needed only
+# for visualization. To obtain reconstructed sensor-space data for further
+# processing, use :meth:`~mne.Evoked.reconstruct_proj` directly. It can also
+# restrict the operation to selected projectors, which is useful when comparing
+# the effects of different artifact corrections. Here, ``evoked_eeg`` contains
+# both ECG and EOG projectors, but we apply and reconstruct using only the EOG
+# projectors:
+
+evoked_eog_reconstructed = evoked_eeg.copy().reconstruct_proj(projs=eog_projs)
+evoked_eog_reconstructed.plot(spatial_colors=True)
+
 # %%
 # Note that here the bias in the EEG and magnetometer channels is reduced by
 # the reconstruction. This suggests that the application of SSP has slightly