This fails the unit handling...

That's tricky: so it seems that the Quantity class being tested is not an ndarray, but tries to behave similarly so it can be used as an array. But it also has an __array__ method as well which converts it to a plain ndarray. I guess the problem is that the presence of __array__ can mean both "I'm not an array, so if you want an array, use this" or "I'm an array duck type, but if you want a real ndarray, use this". I think ideally this whole routine would be something like

return x if isinstance(x, (np.ndarray, ndarray_ducktype)) else np.asarray(x)

But the real question is how to determine whether something is an ndarray_ducktype...

But the real question is how to determine whether something is an ndarray_ducktype ...

That cannot be answered generally. Duck-typing is an operational concept: If it provides the functions I need, it's a duck-type to me. Duck types are never perfect (at the very least they'd fail isinstance() behavior), but if I only need objects that quack like a duck I don't care whether it walks like a duck as well.

The only thing we could do is hard-code a list of known good-enough duck-types. Otherwise we'd have to query the object upfront for all the behavior we want which is practially impossible.
A third theoretical option would be to fully embrace ducktyping: Try with the object as is, if that fails convert the object to a numpy array and retry. However, that can lead to performance issues and more importantly our code is not designed for retries, so we might get into consistency issues if the first try partially modifies some state.

Just writing as I sort this out in my head.

Most of our methods do something like

        x, y = self._process_unit_info([("x", x), ("y", y)], kwargs)
        x = np.ma.ravel(x)  # or alternately np.asarray(x)
        y = np.ma.ravel(y)

by which point we have set the units, and x and y are now type ndarray.

The ax.plot method, on the other hand, needs to unpack x and y into individual arguments so we do:

        x = _check_1d(xy[0])
        y = _check_1d(xy[1])
        self.axes.xaxis.update_units(x)
        self.axes.yaxis.update_units(y)

However, we do not want to lose units yet, and indeed we do not, at least according to test_units.py::test_numpy_facade, which claims we have the ability to transform units at any time with ax.xaxis.set_units. However, note that ax.xaxis.set_units absolutely fails with scatter (see #9713 (comment)) so dragging around the unit information only works for plot currently.

Hence we are in a bit of a Gordian knot of compatibility. We have enshrined pint as something we support, but only for plot apparently. We allow __array__ but currently depend on it having ndim and shape. If an array doesn't have shape then we np.atleast_1d it. I suppose we can do the same thing to arrays without ndim.

This works, but I really feel this is all a mess. Hopefully one of our RSE's can work on this...

@jklymak - indeed, also for astropy unit support works really well for plots, but not for images (cannot give units for an extent). Would be nice if that could also be exchanged. I think for the particular bug that triggered all this, we can also adjust on our side -- it is not something used all the time.

Just to clarify my comment, again while I am thinking about it - ideally plot would do the same thing as the other methods and convert to numpy immediately after the units have been determined. But currently we cannot because plot (and only plot) tries to preserve unit information.

In my mind what the proper way to do this would be to

1. determine units of the object and pass to the axis
2. use the units to covert the unit-full object to numpy into some base unit for the axis (so if the converter allows 'feet', 'miles', and 'meters', convert to floats in meters).
3. Allow the axis to scale its data in the units that the user wants (so if they specify 'feet', add a transform from meters to feet to the displayed data, but keep the underlying data in meters).

This way Matplotlib wouldn't need to drag around the extra units internally, it would just have a transform on each axis.

yes, that makes sense

I think the reason it's currently set up this way with Pint and plot() is because plot() is the only one I've personally needed to adjust axis units--whereas most of my use of other plot units is spatially, where units only matter for the z/image axis (e.g. imshow, contour).

I'm not sure how you can tweak the current conversion interface without breaking backwards compatibility...

I think my proposal would keep the same public user interface but change the internal representation of data passed in w units.

I think it could also be done so that most features work for downstream libraries the way they do now for scatter (not able to change axis units properly) but to have the plot-like full feature downstream libraries would need to add more logic to their locators and formatters.

See #23015 for more discussion of the convert-immediately strategy for dealing with units.

Despite the discussion above, I think this should go in as-is for 3.5.3....