diff --git a/doc/api/next_api_changes/behavior/19368-DS.rst b/doc/api/next_api_changes/behavior/19368-DS.rst
new file mode 100644
index 000000000000..455eba1570f9
--- /dev/null
+++ b/doc/api/next_api_changes/behavior/19368-DS.rst
@@ -0,0 +1,8 @@
+Large ``imshow`` images are now downsampled
+-------------------------------------------
+When showing an image using `~matplotlib.axes.Axes.imshow` that has more than
+:math:`2^{24}` columns or :math:`2^{23}` rows, the image will now be downsampled
+to below this resolution before being resampled for display by the AGG renderer.
+Previously such a large image would be shown incorrectly. To prevent this
+downsampling and the warning it raises, manually downsample your data before
+handing it to `~matplotlib.axes.Axes.imshow`
diff --git a/lib/matplotlib/image.py b/lib/matplotlib/image.py
index 436589f30dd1..8e89a3b065fe 100644
--- a/lib/matplotlib/image.py
+++ b/lib/matplotlib/image.py
@@ -7,6 +7,7 @@
 import os
 import logging
 from pathlib import Path
+import warnings
 
 import numpy as np
 import PIL.PngImagePlugin
@@ -166,7 +167,22 @@ def _resample(
     allocating the output array and fetching the relevant properties from the
     Image object *image_obj*.
     """
-
+    # AGG can only handle coordinates smaller than 24-bit signed integers,
+    # so raise errors if the input data is larger than _image.resample can
+    # handle.
+    msg = ('Data with more than {n} cannot be accurately displayed. '
+           'Downsampling to less than {n} before displaying. '
+           'To remove this warning, manually downsample your data.')
+    if data.shape[1] > 2**23:
+        warnings.warn(msg.format(n='2**23 columns'))
+        step = int(np.ceil(data.shape[1] / 2**23))
+        data = data[:, ::step]
+        transform = Affine2D().scale(step, 1) + transform
+    if data.shape[0] > 2**24:
+        warnings.warn(msg.format(n='2**24 rows'))
+        step = int(np.ceil(data.shape[0] / 2**24))
+        data = data[::step, :]
+        transform = Affine2D().scale(1, step) + transform
     # decide if we need to apply anti-aliasing if the data is upsampled:
     # compare the number of displayed pixels to the number of
     # the data pixels.
diff --git a/lib/matplotlib/tests/test_image.py b/lib/matplotlib/tests/test_image.py
index 719b19057875..396b7838cd36 100644
--- a/lib/matplotlib/tests/test_image.py
+++ b/lib/matplotlib/tests/test_image.py
@@ -1376,3 +1376,43 @@ def test_rgba_antialias():
     # alternating red and blue stripes become purple
     axs[3].imshow(aa, interpolation='antialiased', interpolation_stage='rgba',
                   cmap=cmap, vmin=-1.2, vmax=1.2)
+
+
+# We check for the warning with a draw() in the test, but we also need to
+# filter the warning as it is emitted by the figure test decorator
+@pytest.mark.filterwarnings(r'ignore:Data with more than .* '
+                            'cannot be accurately displayed')
+@pytest.mark.parametrize('origin', ['upper', 'lower'])
+@pytest.mark.parametrize(
+    'dim, size, msg', [['row', 2**23, r'2\*\*23 columns'],
+                       ['col', 2**24, r'2\*\*24 rows']])
+@check_figures_equal(extensions=('png', ))
+def test_large_image(fig_test, fig_ref, dim, size, msg, origin):
+    # Check that Matplotlib downsamples images that are too big for AGG
+    # See issue #19276. Currently the fix only works for png output but not
+    # pdf or svg output.
+    ax_test = fig_test.subplots()
+    ax_ref = fig_ref.subplots()
+
+    array = np.zeros((1, size + 2))
+    array[:, array.size // 2:] = 1
+    if dim == 'col':
+        array = array.T
+    im = ax_test.imshow(array, vmin=0, vmax=1,
+                        aspect='auto', extent=(0, 1, 0, 1),
+                        interpolation='none',
+                        origin=origin)
+
+    with pytest.warns(UserWarning,
+                      match=f'Data with more than {msg} cannot be '
+                      'accurately displayed.'):
+        fig_test.canvas.draw()
+
+    array = np.zeros((1, 2))
+    array[:, 1] = 1
+    if dim == 'col':
+        array = array.T
+    im = ax_ref.imshow(array, vmin=0, vmax=1, aspect='auto',
+                       extent=(0, 1, 0, 1),
+                       interpolation='none',
+                       origin=origin)