As a result, I am unsure how to write a test for the astropy bug here, though I could easily add a test in astropy itself (which is regularily tested against matplotlib dev).

Can you create a fake unit that mimics the astropy behavior? Something like a very thin wrapper over nd.array that triggers the error but is fixed here? Something like #26953

Thank you for your suggestion ! I wasn't able to construct such a mock with a bottom-up approach (build one from the ground up), supposedly because of the sheer complexity of astropy.units.Quantity, so I'm now trying a top-down approach (take astropy and iteratively remove parts that don't play a role in the error). It's taking a while but I'm getting there.

so I'm now trying a top-down approach (take astropy and iteratively remove parts that don't play a role in the error). It's taking a while but I'm getting there.

Is much appreciated since that helps us better understand what's going wonky here.

Ok I think it's now ready for review.

Thanks!

Re: mocking astropy.units.Quantity -- @mhvk can double check if needed.

Apologies for asking this several hours late, but are either of these relevant here?

@rcomer , not sure. I tracked that to #9049 and only pint was mentioned, not astropy.units.

And the second one to #18286 that mentioned unyt and not astropy.units.

@rcomer Indeed these classes will not reproduce the issue we're currently seeing with astropy.units.Quantity. I think I may have used too broad a title for this PR, but it's really about some corner case specific to astropy.units.Quantity. For instance, astropy.units.Quantity intentionally overrides the ndarray.item method so that it doesn't return numerical scalars, but if I remove def item(...) from my reduced mock, the test passes without a patch.

but if I remove def item(...) from my reduced mock, the test passes without a patch.

Can you add a comment that that's the important bit?

Actually I think every part of the mock class is crucial to reproduce the problem, should I really comment all of them ?
In the real-life use case (with astropy.units.Quantity), the error we see is raised through three nested try/except blocks 🙃

In the real-life use case (with astropy.units.Quantity), the error we see is raised through three nested try/except blocks

😬 No, if it's a combination of the things then yeah it doesn't need a specific comment.

I've been trying to understand what's going on here. It seems:

• _make_image(A, ...) is given A as a unitful array
• a_min = A.min() / a_max = A.max() therefore have units
• so the fix in this PR is to just strip the units by calling np.float64

I think the better long term fix here is to add proper support for the image data to have units in the _ImageBase class. This would provide a path to solving other bugs and feature requests around images with units, e.g. #25062, #19476.

So it's unclear to me whether it's worth merging this PR as a bit of a stopgap fix, or whether we should be aiming to add proper support for unitful data in _ImageBase instead. It would be good to hear what @ksunden thinks on this front from a Matplotlib/units point of view?

As a quick proof of concept, the following diff fixes the original issue, also allows mousing over the data to work (that's still broken with this PR), and opens the door to further improvements like labelling the colorbar with the units.

diff --git a/lib/matplotlib/image.py b/lib/matplotlib/image.py
index e326906269..88325a59ed 100644
--- a/lib/matplotlib/image.py
+++ b/lib/matplotlib/image.py
@@ -26,6 +26,7 @@ import matplotlib.colors as mcolors
 from matplotlib.transforms import (
     Affine2D, BboxBase, Bbox, BboxTransform, BboxTransformTo,
     IdentityTransform, TransformedBbox)
+import matplotlib.units as munits
 
 _log = logging.getLogger(__name__)
 
@@ -724,6 +725,13 @@ class _ImageBase(martist.Artist, cm.ScalarMappable):
         ----------
         A : array-like or `PIL.Image.Image`
         """
+        converter = munits.registry.get_converter(A)
+        if converter is not None:
+            self._units = converter.default_units(A, self)
+            A = converter.convert(A, self._units, self)
+        else:
+            self._units = None
+
         if isinstance(A, PIL.Image.Image):
             A = pil_to_array(A)  # Needed e.g. to apply png palette.
         self._A = self._normalize_image_array(A)

It does have the issue that the last arguments to the converter methods are meant to be of type Axis and not _ImageBase though, but in reality I'm not sure that's a huge issue...

Happy to drop this PR in favour of a better fix if needed !

We don't have converter logic for data passed to color mapping. Folks are asked to pass a numpy array to scalar mappables. Maybe its OK to normalize those mappable to so something without units earlier, but I don't know that we have enough infrastructure in place to actually use unit converters here.

but I don't know that we have enough infrastructure in place to actually use unit converters here.

I still think that throwing a convertor on NoNorm (or possibly writing a UnitNorm() if we want to avoid any magic) might do the trick?

I think the complicated part has always been propagating the correct labels back to the colorbar.