diff --git a/doc/changes/devel.rst b/doc/changes/devel.rst
index d5089fd95ee..f6108f65b30 100644
--- a/doc/changes/devel.rst
+++ b/doc/changes/devel.rst
@@ -56,6 +56,7 @@ Bugs
 - Fix bug where ``encoding`` argument was ignored when reading annotations from an EDF file (:gh:`11958` by :newcontrib:`Andrew Gilbert`)
 - Mark tests ``test_adjacency_matches_ft`` and ``test_fetch_uncompressed_file`` as network tests (:gh:`12041` by :newcontrib:`Maksym Balatsko`)
 - Fix bug with :func:`mne.channels.read_ch_adjacency` (:gh:`11608` by :newcontrib:`Ivan Zubarev`)
+- Fix bug where ``epochs.get_data(..., scalings=...)`` would errantly modify the preloaded data (:gh:`12121` by :newcontrib:`Pablo Mainar` and `Eric Larson`_)
 - Fix bugs with saving splits for :class:`~mne.Epochs` (:gh:`11876` by `Dmitrii Altukhov`_)
 - Fix bug with multi-plot 3D rendering where only one plot was updated (:gh:`11896` by `Eric Larson`_)
 - Fix bug where ``verbose`` level was not respected inside parallel jobs (:gh:`12154` by `Eric Larson`_)
@@ -92,6 +93,7 @@ Bugs
 
 API changes
 ~~~~~~~~~~~
+- The default for :meth:`mne.Epochs.get_data` of ``copy=False`` will change to ``copy=True`` in 1.7. Set it explicitly to avoid a warning (:gh:`12121` by :newcontrib:`Pablo Mainar` and `Eric Larson`_)
 - ``mne.preprocessing.apply_maxfilter`` and ``mne maxfilter`` have been deprecated and will be removed in 1.7. Use :func:`mne.preprocessing.maxwell_filter` (see :ref:`this tutorial <tut-artifact-sss>`) in Python or the command-line utility from MEGIN ``maxfilter`` and :func:`mne.bem.fit_sphere_to_headshape` instead (:gh:`11938` by `Eric Larson`_)
 - :func:`mne.io.kit.read_mrk` reading pickled files is deprecated using something like ``np.savetxt(fid, pts, delimiter="\t", newline="\n")`` to save your points instead (:gh:`11937` by `Eric Larson`_)
 - Replace legacy ``inst.pick_channels`` and ``inst.pick_types`` with ``inst.pick`` (where ``inst`` is an instance of :class:`~mne.io.Raw`, :class:`~mne.Epochs`, or :class:`~mne.Evoked`) wherever possible (:gh:`11907` by `Clemens Brunner`_)
diff --git a/doc/changes/names.inc b/doc/changes/names.inc
index 925e38bfd22..3d441e2289f 100644
--- a/doc/changes/names.inc
+++ b/doc/changes/names.inc
@@ -406,6 +406,8 @@
 
 .. _Pablo-Arias: https://github.com/Pablo-Arias
 
+.. _Pablo Mainar: https://github.com/pablomainar
+
 .. _Padma Sundaram: https://www.nmr.mgh.harvard.edu/user/8071
 
 .. _Paul Pasler: https://github.com/ppasler
diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index 51c6c847de6..f71a07ce7d5 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -49,6 +49,7 @@
 
 # %%
 # The data covariance has an interesting structure because of densely packed sensors:
+
 cov = mne.compute_covariance(epochs, tmax=-0.01)
 mne.viz.plot_cov(cov, raw.info)
 
diff --git a/examples/datasets/limo_data.py b/examples/datasets/limo_data.py
index 62fb0322295..4f5291358c3 100644
--- a/examples/datasets/limo_data.py
+++ b/examples/datasets/limo_data.py
@@ -37,7 +37,6 @@
 # License: BSD-3-Clause
 
 # %%
-
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/decoding/decoding_csp_eeg.py b/examples/decoding/decoding_csp_eeg.py
index 896489252f2..9be079d928f 100644
--- a/examples/decoding/decoding_csp_eeg.py
+++ b/examples/decoding/decoding_csp_eeg.py
@@ -78,8 +78,8 @@
 
 # Define a monte-carlo cross-validation generator (reduce variance):
 scores = []
-epochs_data = epochs.get_data()
-epochs_data_train = epochs_train.get_data()
+epochs_data = epochs.get_data(copy=False)
+epochs_data_train = epochs_train.get_data(copy=False)
 cv = ShuffleSplit(10, test_size=0.2, random_state=42)
 cv_split = cv.split(epochs_data_train)
 
diff --git a/examples/decoding/decoding_csp_timefreq.py b/examples/decoding/decoding_csp_timefreq.py
index c1f88588326..cfeaf326ce6 100644
--- a/examples/decoding/decoding_csp_timefreq.py
+++ b/examples/decoding/decoding_csp_timefreq.py
@@ -105,7 +105,7 @@
     epochs.drop_bad()
     y = le.fit_transform(epochs.events[:, 2])
 
-    X = epochs.get_data()
+    X = epochs.get_data(copy=False)
 
     # Save mean scores over folds for each frequency and time window
     freq_scores[freq] = np.mean(
@@ -165,7 +165,7 @@
         w_tmax = w_time + w_size / 2.0
 
         # Crop data into time-window of interest
-        X = epochs.copy().crop(w_tmin, w_tmax).get_data()
+        X = epochs.get_data(tmin=w_tmin, tmax=w_tmax, copy=False)
 
         # Save mean scores over folds for each frequency and time window
         tf_scores[freq, t] = np.mean(
diff --git a/examples/decoding/decoding_time_generalization_conditions.py b/examples/decoding/decoding_time_generalization_conditions.py
index beb69831b8b..a81123c73d7 100644
--- a/examples/decoding/decoding_time_generalization_conditions.py
+++ b/examples/decoding/decoding_time_generalization_conditions.py
@@ -77,12 +77,12 @@
 
 # Fit classifiers on the epochs where the stimulus was presented to the left.
 # Note that the experimental condition y indicates auditory or visual
-time_gen.fit(X=epochs["Left"].get_data(), y=epochs["Left"].events[:, 2] > 2)
+time_gen.fit(X=epochs["Left"].get_data(copy=False), y=epochs["Left"].events[:, 2] > 2)
 
 # %%
 # Score on the epochs where the stimulus was presented to the right.
 scores = time_gen.score(
-    X=epochs["Right"].get_data(), y=epochs["Right"].events[:, 2] > 2
+    X=epochs["Right"].get_data(copy=False), y=epochs["Right"].events[:, 2] > 2
 )
 
 # %%
diff --git a/examples/decoding/decoding_unsupervised_spatial_filter.py b/examples/decoding/decoding_unsupervised_spatial_filter.py
index 07c18813ab8..b27c16e003a 100644
--- a/examples/decoding/decoding_unsupervised_spatial_filter.py
+++ b/examples/decoding/decoding_unsupervised_spatial_filter.py
@@ -58,7 +58,7 @@
     verbose=False,
 )
 
-X = epochs.get_data()
+X = epochs.get_data(copy=False)
 
 ##############################################################################
 # Transform data with PCA computed on the average ie evoked response
diff --git a/examples/decoding/ems_filtering.py b/examples/decoding/ems_filtering.py
index ff6296bf1f1..0273643c61a 100644
--- a/examples/decoding/ems_filtering.py
+++ b/examples/decoding/ems_filtering.py
@@ -64,7 +64,7 @@
 epochs.pick("grad")
 
 # Setup the data to use it a scikit-learn way:
-X = epochs.get_data()  # The MEG data
+X = epochs.get_data(copy=False)  # The MEG data
 y = epochs.events[:, 2]  # The conditions indices
 n_epochs, n_channels, n_times = X.shape
 
diff --git a/examples/decoding/linear_model_patterns.py b/examples/decoding/linear_model_patterns.py
index 4b23e5d1e56..caf53603f97 100644
--- a/examples/decoding/linear_model_patterns.py
+++ b/examples/decoding/linear_model_patterns.py
@@ -60,7 +60,7 @@
 
 # get MEG data
 meg_epochs = epochs.copy().pick(picks="meg", exclude="bads")
-meg_data = meg_epochs.get_data().reshape(len(labels), -1)
+meg_data = meg_epochs.get_data(copy=False).reshape(len(labels), -1)
 
 # %%
 # Decoding in sensor space using a LogisticRegression classifier
diff --git a/examples/decoding/ssd_spatial_filters.py b/examples/decoding/ssd_spatial_filters.py
index a2bdcabf9a1..c3165a7110f 100644
--- a/examples/decoding/ssd_spatial_filters.py
+++ b/examples/decoding/ssd_spatial_filters.py
@@ -146,7 +146,7 @@
         h_trans_bandwidth=1,
     ),
 )
-ssd_epochs.fit(X=epochs.get_data())
+ssd_epochs.fit(X=epochs.get_data(copy=False))
 
 # Plot topographies.
 pattern_epochs = mne.EvokedArray(data=ssd_epochs.patterns_[:4].T, info=ssd_epochs.info)
diff --git a/examples/preprocessing/otp.py b/examples/preprocessing/otp.py
index 7e5e28561fb..afef134c61d 100644
--- a/examples/preprocessing/otp.py
+++ b/examples/preprocessing/otp.py
@@ -13,7 +13,6 @@
 # License: BSD-3-Clause
 
 # %%
-
 import numpy as np
 
 import mne
@@ -70,7 +69,7 @@ def compute_bias(raw):
     sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=None, verbose=False)
     cov = mne.compute_covariance(epochs, tmax=0, method="oas", rank=None, verbose=False)
     idx = epochs.time_as_index(0.036)[0]
-    data = epochs.get_data()[:, :, idx].T
+    data = epochs.get_data(copy=False)[:, :, idx].T
     evoked = mne.EvokedArray(data, epochs.info, tmin=0.0)
     dip = fit_dipole(evoked, cov, sphere, n_jobs=None, verbose=False)[0]
     actual_pos = mne.dipole.get_phantom_dipoles()[0][dipole_number - 1]
diff --git a/examples/stats/sensor_regression.py b/examples/stats/sensor_regression.py
index 4d5b02782f3..28d63360776 100644
--- a/examples/stats/sensor_regression.py
+++ b/examples/stats/sensor_regression.py
@@ -18,10 +18,6 @@
 of the words for which we have EEG activity.
 
 For the general methodology, see e.g. :footcite:`HaukEtAl2006`.
-
-References
-----------
-.. footbibliography::
 """
 # Authors: Tal Linzen <linzen@nyu.edu>
 #          Denis A. Engemann <denis.engemann@gmail.com>
@@ -43,7 +39,7 @@
 epochs = mne.read_epochs(path)
 print(epochs.metadata.head())
 
-##############################################################################
+# %%
 # Psycholinguistically relevant word characteristics are continuous. I.e.,
 # concreteness or imaginability is a graded property. In the metadata,
 # we have concreteness ratings on a 5-point scale. We can show the dependence
@@ -59,7 +55,7 @@
     evokeds, colors=colors, split_legend=True, cmap=(name + " Percentile", "viridis")
 )
 
-##############################################################################
+# %%
 # We observe that there appears to be a monotonic dependence of EEG on
 # concreteness. We can also conduct a continuous analysis: single-trial level
 # regression with concreteness as a continuous (although here, binned)
@@ -72,7 +68,7 @@
         title=cond, ts_args=dict(time_unit="s"), topomap_args=dict(time_unit="s")
     )
 
-##############################################################################
+# %%
 # Because the :func:`~mne.stats.linear_regression` function also estimates
 # p values, we can --
 # after applying FDR correction for multiple comparisons -- also visualise the
@@ -85,3 +81,8 @@
 reject_H0, fdr_pvals = fdr_correction(res["Concreteness"].p_val.data)
 evoked = res["Concreteness"].beta
 evoked.plot_image(mask=reject_H0, time_unit="s")
+
+# %%
+# References
+# ----------
+# .. footbibliography::
diff --git a/mne/_fiff/tests/test_reference.py b/mne/_fiff/tests/test_reference.py
index 4b94b4a8665..3bd540779ed 100644
--- a/mne/_fiff/tests/test_reference.py
+++ b/mne/_fiff/tests/test_reference.py
@@ -620,12 +620,10 @@ def test_add_reference():
     assert_equal(epochs_ref._data.shape[1], epochs._data.shape[1] + 1)
     _check_channel_names(epochs_ref, "Ref")
     ref_idx = epochs_ref.ch_names.index("Ref")
-    ref_data = epochs_ref.get_data()[:, ref_idx, :]
+    ref_data = epochs_ref.get_data(picks=[ref_idx])[:, 0]
     assert_array_equal(ref_data, 0)
     picks_eeg = pick_types(epochs.info, meg=False, eeg=True)
-    assert_array_equal(
-        epochs.get_data()[:, picks_eeg, :], epochs_ref.get_data()[:, picks_eeg, :]
-    )
+    assert_array_equal(epochs.get_data(picks_eeg), epochs_ref.get_data(picks_eeg))
 
     # add two reference channels to epochs
     raw = read_raw_fif(fif_fname, preload=True)
@@ -650,12 +648,10 @@ def test_add_reference():
     ref_idy = epochs_ref.ch_names.index("M2")
     assert_equal(epochs_ref.info["chs"][ref_idx]["ch_name"], "M1")
     assert_equal(epochs_ref.info["chs"][ref_idy]["ch_name"], "M2")
-    ref_data = epochs_ref.get_data()[:, [ref_idx, ref_idy], :]
+    ref_data = epochs_ref.get_data([ref_idx, ref_idy])
     assert_array_equal(ref_data, 0)
     picks_eeg = pick_types(epochs.info, meg=False, eeg=True)
-    assert_array_equal(
-        epochs.get_data()[:, picks_eeg, :], epochs_ref.get_data()[:, picks_eeg, :]
-    )
+    assert_array_equal(epochs.get_data(picks_eeg), epochs_ref.get_data(picks_eeg))
 
     # add reference channel to evoked
     raw = read_raw_fif(fif_fname, preload=True)
@@ -725,7 +721,7 @@ def test_add_reference():
     data = data.get_data()
     epochs = make_fixed_length_epochs(raw).load_data()
     data_2 = epochs.copy().add_reference_channels(["REF"]).pick(picks="eeg")
-    data_2 = data_2.get_data()[0]
+    data_2 = data_2.get_data(copy=False)[0]
     assert_allclose(data, data_2)
     evoked = epochs.average()
     data_3 = evoked.copy().add_reference_channels(["REF"]).pick(picks="eeg")
diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py
index 5fe73244485..a368b7fce0a 100644
--- a/mne/beamformer/_dics.py
+++ b/mne/beamformer/_dics.py
@@ -493,7 +493,7 @@ def apply_dics_epochs(epochs, filters, return_generator=False, verbose=None):
     tmin = epochs.times[0]
 
     sel = _check_channels_spatial_filter(epochs.ch_names, filters)
-    data = epochs.get_data()[:, sel, :]
+    data = epochs.get_data(sel)
 
     stcs = _apply_dics(data=data, filters=filters, info=info, tmin=tmin)
 
diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py
index c096582cc1a..d89fbc35342 100644
--- a/mne/beamformer/_lcmv.py
+++ b/mne/beamformer/_lcmv.py
@@ -402,7 +402,7 @@ def apply_lcmv_epochs(epochs, filters, *, return_generator=False, verbose=None):
     tmin = epochs.times[0]
 
     sel = _check_channels_spatial_filter(epochs.ch_names, filters)
-    data = epochs.get_data()[:, sel, :]
+    data = epochs.get_data(sel)
     stcs = _apply_lcmv(data=data, filters=filters, info=info, tmin=tmin)
 
     if not return_generator:
diff --git a/mne/channels/channels.py b/mne/channels/channels.py
index 725aa884493..85484fa26f1 100644
--- a/mne/channels/channels.py
+++ b/mne/channels/channels.py
@@ -1898,7 +1898,10 @@ def combine_channels(
     ch_idx = list(range(inst.info["nchan"]))
     ch_names = inst.info["ch_names"]
     ch_types = inst.get_channel_types()
-    inst_data = inst.data if isinstance(inst, Evoked) else inst.get_data()
+    kwargs = dict()
+    if isinstance(inst, BaseEpochs):
+        kwargs["copy"] = False
+    inst_data = inst.get_data(**kwargs)
     groups = OrderedDict(deepcopy(groups))
 
     # Convert string values of ``method`` into callables
diff --git a/mne/channels/tests/test_channels.py b/mne/channels/tests/test_channels.py
index 7e27f301048..8e3e482659c 100644
--- a/mne/channels/tests/test_channels.py
+++ b/mne/channels/tests/test_channels.py
@@ -615,7 +615,7 @@ def test_equalize_channels():
     assert raw2.ch_names == ["CH1", "CH2"]
     assert_array_equal(raw2.get_data(), [[1.0], [2.0]])
     assert epochs2.ch_names == ["CH1", "CH2"]
-    assert_array_equal(epochs2.get_data(), [[[3.0], [2.0]]])
+    assert_array_equal(epochs2.get_data(copy=False), [[[3.0], [2.0]]])
     assert cov2.ch_names == ["CH1", "CH2"]
     assert cov2["bads"] == cov["bads"]
     assert ave2.ch_names == ave.ch_names
diff --git a/mne/channels/tests/test_interpolation.py b/mne/channels/tests/test_interpolation.py
index 58cb6a1e669..4f37494e652 100644
--- a/mne/channels/tests/test_interpolation.py
+++ b/mne/channels/tests/test_interpolation.py
@@ -232,10 +232,10 @@ def test_interpolation_meg():
     assert len(raw_meg.info["bads"]) == len(raw_meg.info["bads"])
 
     # MEG -- epochs
-    data1 = epochs_meg.get_data()[:, pick, :].ravel()
+    data1 = epochs_meg.get_data(pick).ravel()
     epochs_meg.info.normalize_proj()
     epochs_meg.interpolate_bads(mode="fast")
-    data2 = epochs_meg.get_data()[:, pick, :].ravel()
+    data2 = epochs_meg.get_data(pick).ravel()
     assert np.corrcoef(data1, data2)[0, 1] > thresh
     assert len(epochs_meg.info["bads"]) == 0
 
diff --git a/mne/decoding/tests/test_base.py b/mne/decoding/tests/test_base.py
index 09dc43c1f8d..885d7ff04f5 100644
--- a/mne/decoding/tests/test_base.py
+++ b/mne/decoding/tests/test_base.py
@@ -317,7 +317,7 @@ def test_get_coef_multiclass_full(n_classes, n_channels, n_times):
     )
     scorer = "roc_auc_ovr_weighted"
     time_gen = GeneralizingEstimator(clf, scorer, verbose=True)
-    X = epochs.get_data()
+    X = epochs.get_data(copy=False)
     y = epochs.events[:, 2]
     n_splits = 3
     cv = StratifiedKFold(n_splits=n_splits)
diff --git a/mne/decoding/tests/test_csp.py b/mne/decoding/tests/test_csp.py
index 788a4040de8..5ad66885392 100644
--- a/mne/decoding/tests/test_csp.py
+++ b/mne/decoding/tests/test_csp.py
@@ -123,7 +123,7 @@ def test_csp():
         preload=True,
         proj=False,
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     n_channels = epochs_data.shape[1]
     y = epochs.events[:, -1]
 
@@ -182,7 +182,7 @@ def test_csp():
         proj=False,
         preload=True,
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     n_channels = epochs_data.shape[1]
 
     n_channels = epochs_data.shape[1]
@@ -256,7 +256,7 @@ def test_regularized_csp():
     epochs = Epochs(
         raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     n_channels = epochs_data.shape[1]
 
     n_components = 3
diff --git a/mne/decoding/tests/test_ems.py b/mne/decoding/tests/test_ems.py
index 6238058658d..44bb0f86135 100644
--- a/mne/decoding/tests/test_ems.py
+++ b/mne/decoding/tests/test_ems.py
@@ -76,7 +76,7 @@ def test_ems():
     raw.close()
 
     # EMS transformer, check that identical to compute_ems
-    X = epochs.get_data()
+    X = epochs.get_data(copy=False)
     y = epochs.events[:, 2]
     X = X / np.std(X)  # X scaled outside cv in compute_ems
     Xt, coefs = list(), list()
diff --git a/mne/decoding/tests/test_transformer.py b/mne/decoding/tests/test_transformer.py
index f1a84c5d41d..cbd586601db 100644
--- a/mne/decoding/tests/test_transformer.py
+++ b/mne/decoding/tests/test_transformer.py
@@ -55,7 +55,7 @@ def test_scaler(info, method):
     epochs = Epochs(
         raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     y = epochs.events[:, -1]
 
     epochs_data_t = epochs_data.transpose([1, 0, 2])
@@ -115,7 +115,7 @@ def test_scaler(info, method):
         picks=np.arange(len(raw.ch_names)),
     )  # non-data chs
     scaler = Scaler(epochs_bad.info, None)
-    pytest.raises(ValueError, scaler.fit, epochs_bad.get_data(), y)
+    pytest.raises(ValueError, scaler.fit, epochs_bad.get_data(copy=False), y)
 
 
 def test_filterestimator():
@@ -129,7 +129,7 @@ def test_filterestimator():
     epochs = Epochs(
         raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
 
     # Add tests for different combinations of l_freq and h_freq
     filt = FilterEstimator(epochs.info, l_freq=40, h_freq=80)
@@ -180,7 +180,7 @@ def test_psdestimator():
     epochs = Epochs(
         raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), preload=True
     )
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     psd = PSDEstimator(2 * np.pi, 0, np.inf)
     y = epochs.events[:, -1]
     X = psd.fit_transform(epochs_data, y)
@@ -244,7 +244,7 @@ def test_unsupervised_spatial_filter():
     pytest.raises(ValueError, UnsupervisedSpatialFilter, KernelRidge(2))
 
     # Test fit
-    X = epochs.get_data()
+    X = epochs.get_data(copy=False)
     n_components = 4
     usf = UnsupervisedSpatialFilter(PCA(n_components))
     usf.fit(X)
diff --git a/mne/epochs.py b/mne/epochs.py
index 50ecf10ee64..b7afada3d1a 100644
--- a/mne/epochs.py
+++ b/mne/epochs.py
@@ -154,7 +154,7 @@ def _save_part(fid, epochs, fmt, n_parts, next_fname, next_idx):
     start_block(fid, FIFF.FIFFB_MNE_EPOCHS)
 
     # write events out after getting data to ensure bad events are dropped
-    data = epochs.get_data()
+    data = epochs.get_data(copy=False)
 
     _check_option("fmt", fmt, ["single", "double"])
 
@@ -1590,6 +1590,7 @@ def _get_data(
         units=None,
         tmin=None,
         tmax=None,
+        copy=False,
         on_empty="warn",
         verbose=None,
     ):
@@ -1612,6 +1613,9 @@ def _get_data(
         """
         from .io.base import _get_ch_factors
 
+        if copy is not None:
+            _validate_type(copy, bool, "copy")
+
         # Handle empty epochs
         self._handle_empty(on_empty, "_get_data")
         # if called with 'out=False', the call came from 'drop_bad()'
@@ -1673,19 +1677,23 @@ def _get_data(
         # handle units param only if we are going to return data (out==True)
         if (units is not None) and out:
             ch_factors = _get_ch_factors(self, units, picks)
+        else:
+            ch_factors = None
 
         if self._bad_dropped:
             if not out:
                 return
             if self.preload:
-                data = data[select]
-                if orig_picks is not None:
-                    data = data[:, picks]
-                if units is not None:
-                    data *= ch_factors[:, np.newaxis]
-                if start != 0 or stop != self.times.size:
-                    data = data[..., start:stop]
-                return data
+                return self._data_sel_copy_scale(
+                    data,
+                    select=select,
+                    orig_picks=orig_picks,
+                    picks=picks,
+                    ch_factors=ch_factors,
+                    start=start,
+                    stop=stop,
+                    copy=copy,
+                )
 
             # we need to load from disk, drop, and return data
             detrend_picks = self._detrend_picks
@@ -1777,16 +1785,63 @@ def _get_data(
                 good_idx, None, copy=False, drop_event_id=False, select_data=False
             )
 
-        if out:
-            if orig_picks is not None:
-                data = data[:, picks]
-            if units is not None:
-                data *= ch_factors[:, np.newaxis]
-            if start != 0 or stop != self.times.size:
-                data = data[..., start:stop]
-            return data
+        if not out:
+            return
+        return self._data_sel_copy_scale(
+            data,
+            select=slice(None),
+            orig_picks=orig_picks,
+            picks=picks,
+            ch_factors=ch_factors,
+            start=start,
+            stop=stop,
+            copy=copy,
+        )
+
+    def _data_sel_copy_scale(
+        self, data, *, select, orig_picks, picks, ch_factors, start, stop, copy
+    ):
+        # data arg starts out as self._data when data is preloaded
+        data_is_self_data = bool(self.preload)
+        logger.debug(f"Data is self data: {data_is_self_data}")
+        # only two types of epoch subselection allowed
+        assert isinstance(select, (slice, np.ndarray)), type(select)
+        if not isinstance(select, slice):
+            logger.debug("  Copying, fancy indexed epochs")
+            data_is_self_data = False  # copy (fancy indexing)
+        elif select != slice(None):
+            logger.debug("  Slicing epochs")
+        if orig_picks is not None:
+            logger.debug("  Copying, fancy indexed picks")
+            assert isinstance(picks, np.ndarray), type(picks)
+            data_is_self_data = False  # copy (fancy indexing)
         else:
-            return None
+            picks = slice(None)
+        if not all(isinstance(x, slice) and x == slice(None) for x in (select, picks)):
+            data = data[select][:, picks]
+        del picks
+        if start != 0 or stop != self.times.size:
+            logger.debug("  Slicing time")
+            data = data[..., start:stop]  # view (slice)
+        if ch_factors is not None:
+            if data_is_self_data:
+                logger.debug("  Copying, scale factors applied")
+                data = data.copy()
+                data_is_self_data = False
+            data *= ch_factors[:, np.newaxis]
+        if not data_is_self_data:
+            return data
+        if copy is None:
+            warn(
+                "The current default of copy=False will change to copy=True in 1.7. "
+                "Set the value of copy explicitly to avoid this warning",
+                FutureWarning,
+            )
+            copy = False
+        if copy:
+            logger.debug("  Copying, copy=True")
+            data = data.copy()
+        return data
 
     @property
     def _detrend_picks(self):
@@ -1797,8 +1852,18 @@ def _detrend_picks(self):
         else:
             return []
 
-    @fill_doc
-    def get_data(self, picks=None, item=None, units=None, tmin=None, tmax=None):
+    @verbose
+    def get_data(
+        self,
+        picks=None,
+        item=None,
+        units=None,
+        tmin=None,
+        tmax=None,
+        *,
+        copy=None,
+        verbose=None,
+    ):
         """Get all epochs as a 3D array.
 
         Parameters
@@ -1823,13 +1888,31 @@ def get_data(self, picks=None, item=None, units=None, tmin=None, tmax=None):
             End time of data to get in seconds.
 
             .. versionadded:: 0.24.0
+        copy : bool
+            Whether to return a copy of the object's data, or (if possible) a view.
+            See :ref:`the NumPy docs <numpy:basics.copies-and-views>` for an
+            explanation. Default is ``False`` in 1.6 but will change to ``True`` in 1.7,
+            set it explicitly to avoid a warning in some cases. A view is only possible
+            when ``item is None``, ``picks is None``, ``units is None``, and data are
+            preloaded.
+
+            .. warning::
+               Using ``copy=False`` and then modifying the returned ``data`` will in
+               turn modify the Epochs object. Use with caution!
+
+            .. versionchanged:: 1.7
+               The default changed from ``False`` to ``True``.
+        %(verbose)s
 
         Returns
         -------
         data : array of shape (n_epochs, n_channels, n_times)
-            A view on epochs data.
+            The epochs data. Will be a copy when ``copy=True`` and will be a view
+            when possible when ``copy=False``.
         """
-        return self._get_data(picks=picks, item=item, units=units, tmin=tmin, tmax=tmax)
+        return self._get_data(
+            picks=picks, item=item, units=units, tmin=tmin, tmax=tmax, copy=copy
+        )
 
     @verbose
     def apply_function(
@@ -2095,7 +2178,7 @@ def save(
             warn("Saving epochs with no data")
             total_size = 0
         else:
-            d = self[0].get_data()
+            d = self[0].get_data(copy=False)
             # this should be guaranteed by subclasses
             assert d.dtype in (">f8", "<f8", ">c16", "<c16")
             total_size = d.nbytes * len(self)
@@ -4201,7 +4284,7 @@ def _concatenate_epochs(
     if with_data:
         offsets = np.cumsum(offsets)
         for start, stop, epochs in zip(offsets[:-1], offsets[1:], epochs_list):
-            this_data = epochs.get_data()
+            this_data = epochs.get_data(copy=False)
             if data is None:
                 data = np.empty(
                     (offsets[-1], len(out.ch_names), len(out.times)),
diff --git a/mne/forward/_field_interpolation.py b/mne/forward/_field_interpolation.py
index 43196bcddbf..272d58e69f3 100644
--- a/mne/forward/_field_interpolation.py
+++ b/mne/forward/_field_interpolation.py
@@ -5,6 +5,7 @@
 # The computations in this code were primarily derived from Matti Hämäläinen's
 # C code.
 
+import inspect
 from copy import deepcopy
 
 import numpy as np
@@ -260,7 +261,10 @@ def _as_meg_type_inst(inst, ch_type="grad", mode="fast"):
     # compute data by multiplying by the 'gain matrix' from
     # original sensors to virtual sensors
     if hasattr(inst, "get_data"):
-        data = inst.get_data()
+        kwargs = dict()
+        if "copy" in inspect.getfullargspec(inst.get_data).kwonlyargs:
+            kwargs["copy"] = False
+        data = inst.get_data(**kwargs)
     else:
         data = inst.data
 
diff --git a/mne/forward/tests/test_field_interpolation.py b/mne/forward/tests/test_field_interpolation.py
index d7d0cb5ac3c..2f2c179e694 100644
--- a/mne/forward/tests/test_field_interpolation.py
+++ b/mne/forward/tests/test_field_interpolation.py
@@ -318,4 +318,4 @@ def test_as_meg_type_evoked():
     virt_epochs.info.normalize_proj()
     virt_epochs = virt_epochs.as_type("mag")
     assert all(ch.endswith("_v") for ch in virt_epochs.info["ch_names"])
-    assert_allclose(virt_epochs.get_data().mean(0), virt_evoked.data)
+    assert_allclose(virt_epochs.get_data(copy=False).mean(0), virt_evoked.data)
diff --git a/mne/io/eeglab/tests/test_eeglab.py b/mne/io/eeglab/tests/test_eeglab.py
index c8519eb335d..ffe67b9a902 100644
--- a/mne/io/eeglab/tests/test_eeglab.py
+++ b/mne/io/eeglab/tests/test_eeglab.py
@@ -355,7 +355,7 @@ def test_io_set_epochs(fnames):
         epochs2 = read_epochs_eeglab(epochs_fname_onefile)
     # one warning for each read_epochs_eeglab because both files have epochs
     # associated with multiple events
-    assert_array_equal(epochs.get_data(), epochs2.get_data())
+    assert_array_equal(epochs.get_data(copy=False), epochs2.get_data(copy=False))
 
 
 @testing.requires_testing_data
diff --git a/mne/io/fieldtrip/tests/test_fieldtrip.py b/mne/io/fieldtrip/tests/test_fieldtrip.py
index 37b5bb2cd41..fc773f85eb8 100644
--- a/mne/io/fieldtrip/tests/test_fieldtrip.py
+++ b/mne/io/fieldtrip/tests/test_fieldtrip.py
@@ -117,8 +117,8 @@ def test_read_epochs(cur_system, version, use_info, monkeypatch):
             if info is None:
                 assert_warning_in_record(NOINFO_WARNING, warn_record)
 
-    mne_data = mne_epoched.get_data()[:, :, :-1]
-    ft_data = epoched_ft.get_data()
+    mne_data = mne_epoched.get_data(copy=False)[:, :, :-1]
+    ft_data = epoched_ft.get_data(copy=False)
 
     check_data(mne_data, ft_data, cur_system)
     check_info_fields(mne_epoched, epoched_ft, use_info)
diff --git a/mne/io/kit/tests/test_kit.py b/mne/io/kit/tests/test_kit.py
index 6e9543ac0bd..244fa41f56b 100644
--- a/mne/io/kit/tests/test_kit.py
+++ b/mne/io/kit/tests/test_kit.py
@@ -286,9 +286,9 @@ def test_epochs():
     raw = read_raw_kit(sqd_path, stim=None)
     events = read_events(events_path)
     raw_epochs = Epochs(raw, events, None, tmin=0, tmax=0.099, baseline=None)
-    data1 = raw_epochs.get_data()
+    data1 = raw_epochs.get_data(copy=False)
     epochs = read_epochs_kit(epochs_path, events_path)
-    data11 = epochs.get_data()
+    data11 = epochs.get_data(copy=False)
     assert_array_equal(data1, data11)
 
 
diff --git a/mne/minimum_norm/time_frequency.py b/mne/minimum_norm/time_frequency.py
index f9f5571ae9b..a8699327ccf 100644
--- a/mne/minimum_norm/time_frequency.py
+++ b/mne/minimum_norm/time_frequency.py
@@ -516,7 +516,7 @@ def _source_induced_power(
                     "single `label=label`, or set `label=None`."
                 )
 
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     K, sel, Vh, vertno, is_free_ori, noise_norm, k_id = _prepare_source_params(
         inst=epochs,
         inverse_operator=inverse_operator,
diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py
index df1952669a5..9d32eabea53 100644
--- a/mne/preprocessing/ica.py
+++ b/mne/preprocessing/ica.py
@@ -818,7 +818,7 @@ def _fit_epochs(self, epochs, picks, decim, verbose):
             )
 
         # this should be a copy (picks a list of int)
-        data = epochs.get_data()[:, picks]
+        data = epochs.get_data(picks=picks)
         # this will be a view
         if decim is not None:
             data = data[:, :, ::decim]
@@ -1036,7 +1036,7 @@ def _transform_epochs(self, epochs, concatenate):
         if not hasattr(self, "mixing_matrix_"):
             raise RuntimeError("No fit available. Please fit ICA.")
         picks = self._get_picks(epochs)
-        data = np.hstack(epochs.get_data()[:, picks])
+        data = np.hstack(epochs.get_data(picks=picks))
         sources = self._transform(data)
         if not concatenate:
             # Put the data back in 3D
@@ -1478,7 +1478,7 @@ def _check_target(self, target, inst, start, stop, reject_by_annotation=False):
         elif isinstance(inst, BaseEpochs):
             if isinstance(target, str):
                 pick = _get_target_ch(inst, target)
-                target = inst.get_data()[:, pick]
+                target = inst.get_data(picks=pick)
 
             if hasattr(target, "ndim"):
                 if target.ndim == 3 and min(target.shape) == 1:
@@ -1710,7 +1710,7 @@ def find_bads_ecg(
                         keep_ecg=False,
                         reject_by_annotation=reject_by_annotation,
                     )
-                ).get_data()
+                ).get_data(copy=False)
 
                 if sources.shape[0] == 0:
                     warn(
@@ -1718,7 +1718,7 @@ def find_bads_ecg(
                         "the input parameters."
                     )
             elif isinstance(inst, BaseEpochs):
-                sources = self.get_sources(inst).get_data()
+                sources = self.get_sources(inst).get_data(copy=False)
             else:
                 raise ValueError(
                     "With `ctps` only Raw and Epochs input is " "supported"
diff --git a/mne/preprocessing/tests/test_csd.py b/mne/preprocessing/tests/test_csd.py
index 1387b6b4465..5ab63b098ff 100644
--- a/mne/preprocessing/tests/test_csd.py
+++ b/mne/preprocessing/tests/test_csd.py
@@ -150,7 +150,7 @@ def test_csd_degenerate(evoked_csd_sphere):
     )
     epochs.drop_bad()
     assert len(epochs) == 1
-    assert_allclose(epochs.get_data()[0], evoked.data)
+    assert_allclose(epochs.get_data(item=[0])[0], evoked.data)
     with pytest.raises(RuntimeError, match="Computing CSD requires.*preload"):
         compute_current_source_density(epochs)
     epochs.load_data()
@@ -159,7 +159,7 @@ def test_csd_degenerate(evoked_csd_sphere):
     evoked = compute_current_source_density(evoked)
     assert_allclose(raw.get_data(), evoked.data)
     epochs = compute_current_source_density(epochs)
-    assert_allclose(epochs.get_data()[0], evoked.data)
+    assert_allclose(epochs.get_data(item=[0])[0], evoked.data)
 
 
 def test_csd_fif():
diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py
index 4f623e202ff..62aab79a3c2 100644
--- a/mne/preprocessing/tests/test_ica.py
+++ b/mne/preprocessing/tests/test_ica.py
@@ -108,7 +108,7 @@ def test_ica_full_data_recovery(method):
     evoked = epochs.average()
     n_channels = 5
     data = raw._data[:n_channels].copy()
-    data_epochs = epochs.get_data()
+    data_epochs = epochs.get_data(copy=True)
     data_evoked = evoked.data
     raw.set_annotations(Annotations([0.5], [0.5], ["BAD"]))
     methods = [method]
@@ -137,7 +137,7 @@ def test_ica_full_data_recovery(method):
                 ica.fit(epochs, picks=picks)
             _assert_ica_attributes(ica, epochs.get_data(picks))
             epochs2 = ica.apply(epochs.copy(), **kwargs)
-            data2 = epochs2.get_data()[:, :n_channels]
+            data2 = epochs2.get_data(picks=slice(0, n_channels))
             if ok:
                 assert_allclose(
                     data_epochs[:, :n_channels], data2, rtol=1e-10, atol=1e-15
@@ -544,13 +544,13 @@ def test_ica_core(method, n_components, noise_cov, n_pca_components, browser_bac
     ica = ICA(noise_cov=noise_cov, n_components=n_components, method=method)
     with _record_warnings():  # sometimes warns
         ica.fit(epochs)
-    _assert_ica_attributes(ica, epochs.get_data(), limits=(0.2, 20))
-    data = epochs.get_data()[:, 0, :]
+    _assert_ica_attributes(ica, epochs.get_data(copy=False), limits=(0.2, 20))
+    data = epochs.get_data(picks=[0])[:, 0]
     n_samples = np.prod(data.shape)
     assert_equal(ica.n_samples_, n_samples)
     print(ica)  # to test repr
 
-    sources = ica.get_sources(epochs).get_data()
+    sources = ica.get_sources(epochs).get_data(copy=False)
     assert sources.shape[1] == ica.n_components_
 
     with pytest.raises(ValueError, match="target do not have the same nu"):
@@ -878,7 +878,7 @@ def f(x, y):
     evoked = epochs.average()
     evoked_data = evoked.data.copy()
     raw_data = raw[:][0].copy()
-    epochs_data = epochs.get_data().copy()
+    epochs_data = epochs.get_data(copy=True)
 
     with pytest.warns(RuntimeWarning, match="longer"):
         idx, scores = ica.find_bads_ecg(
@@ -930,7 +930,7 @@ def f(x, y):
     assert_equal(len(scores), ica.n_components_)
 
     assert_array_equal(raw_data, raw[:][0])
-    assert_array_equal(epochs_data, epochs.get_data())
+    assert_array_equal(epochs_data, epochs.get_data(copy=False))
     assert_array_equal(evoked_data, evoked.data)
 
     # check score funcs
@@ -984,7 +984,7 @@ def f(x, y):
     assert ica_epochs.events.shape == epochs.events.shape
     ica_chans = [ch for ch in ica_epochs.ch_names if "ICA" in ch]
     assert ica.n_components_ == len(ica_chans)
-    assert ica.n_components_ == ica_epochs.get_data().shape[1]
+    assert ica.n_components_ == ica_epochs.get_data(copy=False).shape[1]
     assert ica_epochs._raw is None
     assert ica_epochs.preload is True
 
@@ -1116,7 +1116,7 @@ def test_ica_cov(method, cov, tmp_path, short_raw_epochs):
     with _record_warnings():  # ICA does not converge
         ica.fit(raw, picks=np.arange(10))
     _assert_ica_attributes(ica)
-    sources = ica.get_sources(epochs).get_data()
+    sources = ica.get_sources(epochs).get_data(copy=False)
     assert ica.mixing_matrix_.shape == (2, 2)
     assert ica.unmixing_matrix_.shape == (2, 2)
     assert ica.pca_components_.shape == (10, 10)
diff --git a/mne/preprocessing/tests/test_regress.py b/mne/preprocessing/tests/test_regress.py
index ab2a640343b..860481207a4 100644
--- a/mne/preprocessing/tests/test_regress.py
+++ b/mne/preprocessing/tests/test_regress.py
@@ -31,7 +31,7 @@ def test_regress_artifact():
     orig_norm = np.linalg.norm(orig_data)
     epochs_clean, betas = regress_artifact(epochs)
     regress_artifact(epochs, betas=betas, copy=False)  # inplace, and w/betas
-    assert_allclose(epochs_clean.get_data(), epochs.get_data())
+    assert_allclose(epochs_clean.get_data(copy=False), epochs.get_data(copy=False))
     clean_data = epochs_clean.get_data("eeg")
     clean_norm = np.linalg.norm(clean_data)
     assert orig_norm / 2 > clean_norm > orig_norm / 10
diff --git a/mne/preprocessing/tests/test_stim.py b/mne/preprocessing/tests/test_stim.py
index a639ad5e5d1..e4934488a45 100644
--- a/mne/preprocessing/tests/test_stim.py
+++ b/mne/preprocessing/tests/test_stim.py
@@ -41,18 +41,18 @@ def test_fix_stim_artifact():
     epochs = fix_stim_artifact(
         epochs, tmin=tmin, tmax=tmax, mode="linear", picks=("eeg", "eog")
     )
-    data = epochs.copy().pick(("eeg", "eog")).get_data()[:, :, tmin_samp:tmax_samp]
+    data = epochs.get_data(("eeg", "eog"))[:, :, tmin_samp:tmax_samp]
     diff_data0 = np.diff(data[0][0])
     diff_data0 -= np.mean(diff_data0)
     assert_array_almost_equal(diff_data0, np.zeros(len(diff_data0)))
 
-    data = epochs.copy().pick(("meg")).get_data()[:, :, tmin_samp:tmax_samp]
+    data = epochs.get_data("meg")[:, :, tmin_samp:tmax_samp]
     diff_data0 = np.diff(data[0][0])
     diff_data0 -= np.mean(diff_data0)
     assert np.all(diff_data0 != 0)
 
     epochs = fix_stim_artifact(epochs, tmin=tmin, tmax=tmax, mode="window")
-    data_from_epochs_fix = epochs.get_data()[:, :, tmin_samp:tmax_samp]
+    data_from_epochs_fix = epochs.get_data(copy=False)[:, :, tmin_samp:tmax_samp]
     assert not np.all(data_from_epochs_fix != 0)
 
     # use window before stimulus in raw
@@ -99,7 +99,7 @@ def test_fix_stim_artifact():
     e_start = int(np.ceil(epochs.info["sfreq"] * epochs.tmin))
     tmin_samp = int(-0.035 * epochs.info["sfreq"]) - e_start
     tmax_samp = int(-0.015 * epochs.info["sfreq"]) - e_start
-    data_from_raw_fix = epochs.get_data()[:, :, tmin_samp:tmax_samp]
+    data_from_raw_fix = epochs.get_data(copy=False)[:, :, tmin_samp:tmax_samp]
     assert np.all(data_from_raw_fix) == 0.0
 
     # use window after stimulus
diff --git a/mne/preprocessing/tests/test_xdawn.py b/mne/preprocessing/tests/test_xdawn.py
index 047a35d75dd..59822799853 100644
--- a/mne/preprocessing/tests/test_xdawn.py
+++ b/mne/preprocessing/tests/test_xdawn.py
@@ -59,8 +59,8 @@ def test_xdawn_picks():
     xd.fit(epochs)
     epochs_out = xd.apply(epochs)["1"]
     assert epochs_out.info["ch_names"] == epochs.ch_names
-    assert not (epochs_out.get_data()[:, 0] != data[:, 0]).any()
-    assert_array_equal(epochs_out.get_data()[:, 1], data[:, 1])
+    assert not (epochs_out.get_data([0])[:, 0] != data[:, 0]).any()
+    assert_array_equal(epochs_out.get_data([1])[:, 0], data[:, 1])
 
 
 def test_xdawn_fit():
@@ -375,7 +375,10 @@ def test_xdawn_decoding_performance():
     )
 
     cv = KFold(n_splits=3, shuffle=False)
-    for pipe, X in ((xdawn_pipe, epochs), (xdawn_trans_pipe, epochs.get_data())):
+    for pipe, X in (
+        (xdawn_pipe, epochs),
+        (xdawn_trans_pipe, epochs.get_data(copy=False)),
+    ):
         predictions = np.empty_like(y, dtype=float)
         for train, test in cv.split(X, y):
             pipe.fit(X[train], y[train])
diff --git a/mne/preprocessing/xdawn.py b/mne/preprocessing/xdawn.py
index ab9684cd07d..aed801068ca 100644
--- a/mne/preprocessing/xdawn.py
+++ b/mne/preprocessing/xdawn.py
@@ -451,7 +451,7 @@ def fit(self, epochs, y=None):
             raise ValueError("epochs must be an Epochs object.")
         picks = _pick_data_channels(epochs.info)
         use_info = pick_info(epochs.info, picks)
-        X = epochs.get_data()[:, picks, :]
+        X = epochs.get_data(picks)
         y = epochs.events[:, 2] if y is None else y
         self.event_id_ = epochs.event_id
 
@@ -525,7 +525,7 @@ def transform(self, inst):
             Spatially filtered signals.
         """  # noqa: E501
         if isinstance(inst, BaseEpochs):
-            X = inst.get_data()
+            X = inst.get_data(copy=False)
         elif isinstance(inst, Evoked):
             X = inst.data
         elif isinstance(inst, np.ndarray):
@@ -636,7 +636,7 @@ def _apply_epochs(self, epochs, include, exclude, event_id, picks):
 
         # special case where epochs come picked but fit was 'unpicked'.
         epochs_dict = dict()
-        data = np.hstack(epochs.get_data()[:, picks])
+        data = np.hstack(epochs.get_data(picks))
 
         for eid in event_id:
             data_r = self._pick_sources(data, include, exclude, eid)
diff --git a/mne/rank.py b/mne/rank.py
index 0cab1f3a563..34284db50de 100644
--- a/mne/rank.py
+++ b/mne/rank.py
@@ -442,8 +442,7 @@ def compute_rank(
                 if isinstance(inst, BaseRaw):
                     data = inst.get_data(picks, reject_by_annotation="omit")
                 else:  # isinstance(inst, BaseEpochs):
-                    data = inst.get_data()[:, picks, :]
-                    data = np.concatenate(data, axis=1)
+                    data = np.concatenate(inst.get_data(picks), axis=1)
                 if proj:
                     data = np.dot(proj_op, data)
                 this_rank = _estimate_rank_meeg_signals(
diff --git a/mne/stats/regression.py b/mne/stats/regression.py
index e005832824b..5240f3c61cb 100644
--- a/mne/stats/regression.py
+++ b/mne/stats/regression.py
@@ -76,7 +76,7 @@ def linear_regression(inst, design_matrix, names=None):
         if [inst.ch_names[p] for p in picks] != inst.ch_names:
             warn("Fitting linear model to non-data or bad channels. " "Check picking")
         msg = "Fitting linear model to epochs"
-        data = inst.get_data()
+        data = inst.get_data(copy=False)
         out = EvokedArray(np.zeros(data.shape[1:]), inst.info, inst.tmin)
     elif isgenerator(inst):
         msg = "Fitting linear model to source estimates (generator input)"
diff --git a/mne/stats/tests/test_regression.py b/mne/stats/tests/test_regression.py
index 190e3ceff87..d36bd75f65b 100644
--- a/mne/stats/tests/test_regression.py
+++ b/mne/stats/tests/test_regression.py
@@ -71,6 +71,12 @@ def test_regression():
         for v1, v2 in zip(lm1[k], lm2[k]):
             assert_array_equal(v1.data, v2.data)
 
+    # Smoke test for fitting on epochs
+    epochs.load_data()
+    with pytest.warns(RuntimeWarning, match="non-data"):
+        linear_regression(epochs, design_matrix)
+    linear_regression(epochs.copy().pick("eeg"), design_matrix)
+
 
 @testing.requires_testing_data
 def test_continuous_regression_no_overlap():
diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py
index 423fe556365..e5ac3892ca8 100644
--- a/mne/tests/test_epochs.py
+++ b/mne/tests/test_epochs.py
@@ -85,6 +85,13 @@
 rng = np.random.RandomState(42)
 
 
+pytestmark = [
+    pytest.mark.filterwarnings(
+        "ignore:The current default of copy=False will change to copy=.*:FutureWarning",
+    ),
+]
+
+
 def _create_epochs_with_annotations():
     """Create test dataset of Epochs with Annotations."""
     # set up a test dataset
@@ -276,7 +283,7 @@ def _get_data(preload=False):
 flat = dict(grad=1e-15, mag=1e-15)
 
 
-def test_get_data():
+def test_get_data_copy():
     """Test the .get_data() method."""
     raw, events, picks = _get_data()
     event_id = {"a/1": 1, "a/2": 2, "b/1": 3, "b/2": 4}
@@ -321,6 +328,25 @@ def test_get_data():
     with pytest.raises(TypeError, match="tmax .* float, None"):
         epochs.get_data(tmin=1, tmax=np.ones(5))
 
+    # Test copy
+    data = epochs.get_data(copy=True)
+    assert not np.shares_memory(data, epochs._data)
+
+    with pytest.warns(FutureWarning, match="The current default of copy=False will"):
+        data = epochs.get_data(verbose="debug")
+    assert np.shares_memory(data, epochs._data)
+    assert data is epochs._data
+    data_orig = data.copy()
+    # picks, item, and units must be None
+    data = epochs.get_data(copy=False, picks=[1])
+    assert not np.shares_memory(data, epochs._data)
+    data = epochs.get_data(copy=False, item=[0])
+    assert not np.shares_memory(data, epochs._data)
+    data = epochs.get_data(copy=False, units=dict(eeg="uV"))
+    assert not np.shares_memory(data, epochs._data)
+    # Make sure we didn't mess up our values
+    assert_allclose(data_orig, epochs._data)
+
 
 def test_hierarchical():
     """Test hierarchical access."""
@@ -1033,7 +1059,7 @@ def test_epochs_baseline_basic(preload, tmp_path):
     epochs = mne.Epochs(raw, events, None, 0, 1e-3, baseline=None, preload=preload)
     epochs.drop_bad()
     epochs_nobl = epochs.copy()
-    epochs_data = epochs.get_data()
+    epochs_data = epochs.get_data(copy=False)
     assert epochs_data.shape == (1, 2, 2)
     expected = data.copy()
     assert_array_equal(epochs_data[0], expected)
diff --git a/mne/tests/test_evoked.py b/mne/tests/test_evoked.py
index 9820dcdb5e3..fe3f41f141c 100644
--- a/mne/tests/test_evoked.py
+++ b/mne/tests/test_evoked.py
@@ -135,7 +135,7 @@ def test_decim():
         expected_times = epochs.times[offset::decim]
         assert_allclose(ev_decim.times, expected_times)
         assert_allclose(ev_ep_decim.times, expected_times)
-        expected_data = epochs.get_data()[:, :, offset::decim].mean(axis=0)
+        expected_data = epochs.get_data(copy=False)[:, :, offset::decim].mean(axis=0)
         assert_allclose(ev_decim.data, expected_data)
         assert_allclose(ev_ep_decim.data, expected_data)
         assert_equal(ev_decim.info["sfreq"], sfreq_new)
@@ -911,7 +911,7 @@ def test_hilbert():
     raw_hilb = raw.apply_hilbert()
     epochs_hilb = epochs.apply_hilbert()
     evoked_hilb = evoked.copy().apply_hilbert()
-    evoked_hilb_2_data = epochs_hilb.get_data().mean(0)
+    evoked_hilb_2_data = epochs_hilb.get_data(copy=False).mean(0)
     assert_allclose(evoked_hilb.data, evoked_hilb_2_data)
     # This one is only approximate because of edge artifacts
     evoked_hilb_3 = Epochs(raw_hilb, events).average()
diff --git a/mne/tests/test_filter.py b/mne/tests/test_filter.py
index 552489f45d1..f2b5ec1b2e7 100644
--- a/mne/tests/test_filter.py
+++ b/mne/tests/test_filter.py
@@ -487,7 +487,7 @@ def test_resample_below_1_sample():
     )
     epochs.resample(1)
     assert len(epochs.times) == 1
-    assert epochs.get_data().shape[2] == 1
+    assert epochs.get_data(copy=False).shape[2] == 1
 
 
 @pytest.mark.slowtest
diff --git a/mne/tests/test_rank.py b/mne/tests/test_rank.py
index f88dd68e282..bde640e276c 100644
--- a/mne/tests/test_rank.py
+++ b/mne/tests/test_rank.py
@@ -5,7 +5,14 @@
 import pytest
 from numpy.testing import assert_array_equal
 
-from mne import compute_raw_covariance, pick_info, pick_types, read_cov, read_evokeds
+from mne import (
+    compute_raw_covariance,
+    make_fixed_length_epochs,
+    pick_info,
+    pick_types,
+    read_cov,
+    read_evokeds,
+)
 from mne._fiff.pick import _picks_by_type
 from mne._fiff.proj import _has_eeg_average_ref_proj
 from mne.cov import prepare_noise_cov
@@ -190,6 +197,19 @@ def test_cov_rank_estimation(rank_method, proj, meg):
             assert rank[ch_type] == expected_rank
 
 
+@pytest.mark.parametrize(
+    "rank_method, proj", [("info", True), ("info", False), (None, True), (None, False)]
+)
+def test_rank_epochs(rank_method, proj):
+    """Test that raw and epochs give the same results in a simple case."""
+    # And a smoke test for epochs
+    raw = read_raw_fif(raw_fname, preload=True)
+    epochs = make_fixed_length_epochs(raw, preload=True, proj=False)
+    rank_raw = compute_rank(raw, rank_method, proj=proj)
+    rank_epochs = compute_rank(epochs, rank_method, proj=proj)
+    assert rank_raw == rank_epochs
+
+
 @pytest.mark.slowtest  # ~3 s apiece on Azure means overall it's slow
 @testing.requires_testing_data
 @pytest.mark.parametrize("fname, rank_orig", ((hp_fif_fname, 120), (mf_fif_fname, 67)))
diff --git a/mne/time_frequency/csd.py b/mne/time_frequency/csd.py
index 78ff9d95bcc..e3499f45e2e 100644
--- a/mne/time_frequency/csd.py
+++ b/mne/time_frequency/csd.py
@@ -717,7 +717,7 @@ def csd_fourier(
     """
     epochs, projs = _prepare_csd(epochs, tmin, tmax, picks, projs)
     return csd_array_fourier(
-        epochs.get_data(),
+        epochs.get_data(copy=False),
         sfreq=epochs.info["sfreq"],
         t0=epochs.tmin,
         fmin=fmin,
@@ -900,7 +900,7 @@ def csd_multitaper(
     """
     epochs, projs = _prepare_csd(epochs, tmin, tmax, picks, projs)
     return csd_array_multitaper(
-        epochs.get_data(),
+        epochs.get_data(copy=False),
         sfreq=epochs.info["sfreq"],
         t0=epochs.tmin,
         fmin=fmin,
@@ -1109,7 +1109,7 @@ def csd_morlet(
     """
     epochs, projs = _prepare_csd(epochs, tmin, tmax, picks, projs)
     return csd_array_morlet(
-        epochs.get_data(),
+        epochs.get_data(copy=False),
         sfreq=epochs.info["sfreq"],
         frequencies=frequencies,
         t0=epochs.tmin,
diff --git a/mne/time_frequency/tests/test_csd.py b/mne/time_frequency/tests/test_csd.py
index 6c306188699..3763e6bdbb4 100644
--- a/mne/time_frequency/tests/test_csd.py
+++ b/mne/time_frequency/tests/test_csd.py
@@ -454,7 +454,7 @@ def test_csd_fourier():
     for (tmin, tmax), as_array in parameters:
         if as_array:
             csd = csd_array_fourier(
-                epochs.get_data(),
+                epochs.get_data(copy=False),
                 sfreq,
                 epochs.tmin,
                 fmin=9,
@@ -510,7 +510,7 @@ def test_csd_multitaper():
     for (tmin, tmax), as_array, adaptive in parameters:
         if as_array:
             csd = csd_array_multitaper(
-                epochs.get_data(),
+                epochs.get_data(copy=False),
                 sfreq,
                 epochs.tmin,
                 adaptive=adaptive,
@@ -578,7 +578,7 @@ def test_csd_morlet():
     for (tmin, tmax), as_array in parameters:
         if as_array:
             csd = csd_array_morlet(
-                epochs.get_data(),
+                epochs.get_data(copy=False),
                 sfreq,
                 freqs,
                 t0=epochs.tmin,
diff --git a/mne/time_frequency/tfr.py b/mne/time_frequency/tfr.py
index 09745e9a1b1..bd3a02865c1 100644
--- a/mne/time_frequency/tfr.py
+++ b/mne/time_frequency/tfr.py
@@ -2930,7 +2930,7 @@ def _get_data(inst, return_itc):
     if not isinstance(inst, (BaseEpochs, Evoked)):
         raise TypeError("inst must be Epochs or Evoked")
     if isinstance(inst, BaseEpochs):
-        data = inst.get_data()
+        data = inst.get_data(copy=False)
     else:
         if return_itc:
             raise ValueError("return_itc must be False for evoked data")
diff --git a/mne/utils/tests/test_mixin.py b/mne/utils/tests/test_mixin.py
index 32c7abe6a32..aa13f705f14 100644
--- a/mne/utils/tests/test_mixin.py
+++ b/mne/utils/tests/test_mixin.py
@@ -29,7 +29,11 @@ def test_decimate():
         epo_3=mne.make_fixed_length_epochs(raw, preload=False).decimate(2).decimate(3),
     )
     for key, other in others.items():
-        assert_allclose(epo.get_data(), other.get_data(), err_msg=key)
+        assert_allclose(
+            epo.get_data(copy=False),
+            other.get_data(copy=False),
+            err_msg=key,
+        )
         assert_allclose(epo.times, other.times, err_msg=key)
     evo = epo.average()
     epo_full = mne.make_fixed_length_epochs(raw, preload=True)
diff --git a/mne/viz/_figure.py b/mne/viz/_figure.py
index 7bd58f4ee9e..7f958657876 100644
--- a/mne/viz/_figure.py
+++ b/mne/viz/_figure.py
@@ -327,7 +327,9 @@ def _load_data(self, start=None, stop=None):
             )
             ix_stop = ix_start + self.mne.n_epochs
             item = slice(ix_start, ix_stop)
-            data = np.concatenate(self.mne.inst.get_data(item=item), axis=-1)
+            data = np.concatenate(
+                self.mne.inst.get_data(item=item, copy=False), axis=-1
+            )
             times = np.arange(start, stop) / self.mne.info["sfreq"]
             return data, times
 
@@ -554,7 +556,7 @@ def _create_epoch_histogram(self):
         """Create peak-to-peak histogram of channel amplitudes."""
         epochs = self.mne.inst
         data = OrderedDict()
-        ptp = np.ptp(epochs.get_data(), axis=2)
+        ptp = np.ptp(epochs.get_data(copy=False), axis=2)
         for ch_type in ("eeg", "mag", "grad"):
             if ch_type in epochs:
                 data[ch_type] = ptp.T[self.mne.ch_types == ch_type].ravel()
diff --git a/mne/viz/ica.py b/mne/viz/ica.py
index 24a9af42b12..12002b1cd56 100644
--- a/mne/viz/ica.py
+++ b/mne/viz/ica.py
@@ -226,7 +226,7 @@ def _plot_ica_properties(
 
     # image and erp
     # we create a new epoch with dropped rows
-    epoch_data = epochs_src.get_data()
+    epoch_data = epochs_src.get_data(copy=False)
     epoch_data = np.insert(
         arr=epoch_data,
         obj=(dropped_indices - np.arange(len(dropped_indices))).astype(int),
@@ -744,7 +744,7 @@ def _prepare_data_ica_properties(inst, ica, reject_by_annotation=True, reject="a
         epochs_src = ica.get_sources(inst)
         dropped_indices = []
         kind = "Epochs"
-    return kind, dropped_indices, epochs_src, epochs_src.get_data()
+    return kind, dropped_indices, epochs_src, epochs_src.get_data(copy=False)
 
 
 def _plot_ica_sources_evoked(evoked, picks, exclude, title, show, ica, labels=None):
diff --git a/mne/viz/tests/test_epochs.py b/mne/viz/tests/test_epochs.py
index 77f45ed1598..d3dd90d224d 100644
--- a/mne/viz/tests/test_epochs.py
+++ b/mne/viz/tests/test_epochs.py
@@ -415,7 +415,7 @@ def test_plot_psd_epochs(epochs):
     fig = spectrum.plot_topomap(bands=[(20, "20 Hz"), (15, 25, "15-25 Hz")])
     # test with a flat channel
     err_str = "for channel %s" % epochs.ch_names[2]
-    epochs.get_data()[0, 2, :] = 0
+    epochs.get_data(copy=False)[0, 2, :] = 0
     for dB in [True, False]:
         with pytest.warns(UserWarning, match=err_str):
             epochs.compute_psd().plot(dB=dB)
diff --git a/mne/viz/topo.py b/mne/viz/topo.py
index 8e363d117e9..3cae18c3ed7 100644
--- a/mne/viz/topo.py
+++ b/mne/viz/topo.py
@@ -1242,7 +1242,7 @@ def plot_topo_image_epochs(
     scale_coeffs = [scalings.get(ch_type, 1) for ch_type in ch_types]
     # scale the data
     epochs._data *= np.array(scale_coeffs)[:, np.newaxis]
-    data = epochs.get_data()
+    data = epochs.get_data(copy=False)
     # get vlims for each channel type
     vlim_dict = dict()
     for ch_type in set(ch_types):
diff --git a/tutorials/epochs/10_epochs_overview.py b/tutorials/epochs/10_epochs_overview.py
index 8ec3af8b065..7778110b6a5 100644
--- a/tutorials/epochs/10_epochs_overview.py
+++ b/tutorials/epochs/10_epochs_overview.py
@@ -312,7 +312,9 @@
 shorter_epochs = epochs.copy().crop(tmin=-0.1, tmax=0.1, include_tmax=True)
 
 for name, obj in dict(Original=epochs, Cropped=shorter_epochs).items():
-    print("{} epochs has {} time samples".format(name, obj.get_data().shape[-1]))
+    print(
+        "{} epochs has {} time samples".format(name, obj.get_data(copy=False).shape[-1])
+    )
 
 # %%
 # Cropping removed part of the baseline. When printing the
diff --git a/tutorials/epochs/15_baseline_regression.py b/tutorials/epochs/15_baseline_regression.py
index 7d62a3cae1b..4f3c3456760 100644
--- a/tutorials/epochs/15_baseline_regression.py
+++ b/tutorials/epochs/15_baseline_regression.py
@@ -142,7 +142,7 @@
     epochs.copy()
     .crop(*baseline)
     .pick([ch])
-    .get_data()  # convert to NumPy array
+    .get_data(copy=False)  # convert to NumPy array
     .mean(axis=-1)  # average across timepoints
     .squeeze()  # only 1 channel, so remove singleton dimension
 )
diff --git a/tutorials/epochs/60_make_fixed_length_epochs.py b/tutorials/epochs/60_make_fixed_length_epochs.py
index a311842312a..108435c092e 100644
--- a/tutorials/epochs/60_make_fixed_length_epochs.py
+++ b/tutorials/epochs/60_make_fixed_length_epochs.py
@@ -94,7 +94,7 @@
 # (for more information on filtering, please see :ref:`tut-filter-resample`).
 
 epochs.load_data().filter(l_freq=8, h_freq=12)
-alpha_data = epochs.get_data()
+alpha_data = epochs.get_data(copy=False)
 
 # %%
 # If desired, separate correlation matrices for each epoch can be obtained.
diff --git a/tutorials/machine-learning/50_decoding.py b/tutorials/machine-learning/50_decoding.py
index fe2addc87f3..a15cff55695 100644
--- a/tutorials/machine-learning/50_decoding.py
+++ b/tutorials/machine-learning/50_decoding.py
@@ -82,7 +82,7 @@
 epochs.pick(picks="meg", exclude="bads")  # remove stim and EOG
 del raw
 
-X = epochs.get_data()  # MEG signals: n_epochs, n_meg_channels, n_times
+X = epochs.get_data(copy=False)  # MEG signals: n_epochs, n_meg_channels, n_times
 y = epochs.events[:, 2]  # target: auditory left vs visual left
 
 # %%
diff --git a/tutorials/stats-sensor-space/20_erp_stats.py b/tutorials/stats-sensor-space/20_erp_stats.py
index 08b5a583ffc..cba1d3bbf0f 100644
--- a/tutorials/stats-sensor-space/20_erp_stats.py
+++ b/tutorials/stats-sensor-space/20_erp_stats.py
@@ -90,8 +90,8 @@
 # In this case, inference is done over items. In the same manner, we could
 # also conduct the test over, e.g., subjects.
 X = [
-    long_words.get_data().transpose(0, 2, 1),
-    short_words.get_data().transpose(0, 2, 1),
+    long_words.get_data(copy=False).transpose(0, 2, 1),
+    short_words.get_data(copy=False).transpose(0, 2, 1),
 ]
 tfce = dict(start=0.4, step=0.4)  # ideally start and step would be smaller
 
diff --git a/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py b/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py
index 7c3e272a933..dda2ab29255 100644
--- a/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py
+++ b/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py
@@ -85,7 +85,7 @@
 # Obtain the data as a 3D matrix and transpose it such that
 # the dimensions are as expected for the cluster permutation test:
 # n_epochs × n_times × n_channels
-X = [epochs[event_name].get_data() for event_name in event_id]
+X = [epochs[event_name].get_data(copy=False) for event_name in event_id]
 X = [np.transpose(x, (0, 2, 1)) for x in X]