Probably these changes should be repeated in errorbar3d, but likely after #18436 since it's currently even more divergent. Though I'm also thinking maybe we should refactor some of these two together.

Let's merge the two PRs as separate improvements and figure out how to reuse the logic in both cases later? I actually already have a followup patch for the 2D case that factors together the x and y handling, probably I can make it handle z too.

This breaks mixed unit error bars like:

x = np.arange(10)
y = [datetime(2020, 5, i * 2 + 1) for i in x]
plt.errorbar(x, y, timedelta(days=0.5))

This breaks mixed unit error bars like:

Ooops, good catch!

Good catch, now fixed by using object arrays whenever the inputs are not already arrays (test included).

I think this has broken astropy (astropy/astropy#11479). Any advise would be greatly appreciated. Thanks!

I'll have a look.

Looks like the following fixes the problem:

diff --git i/lib/matplotlib/axes/_axes.py w/lib/matplotlib/axes/_axes.py
index e906dba6c8..e9e8c3ef25 100644
--- i/lib/matplotlib/axes/_axes.py
+++ w/lib/matplotlib/axes/_axes.py
@@ -3427,13 +3427,17 @@ class Axes(_AxesBase):
                 the note in the main docstring about this parameter's name.
             """
             try:
-                low, high = np.broadcast_to(err, (2, len(data)))
+                np.broadcast_to(err, (2, len(data)))
             except ValueError:
                 raise ValueError(
                     f"'{name}err' (shape: {np.shape(err)}) must be a scalar "
                     f"or a 1D or (2, n) array-like whose shape matches "
                     f"'{name}' (shape: {np.shape(data)})") from None
-            return data - low * ~lolims, data + high * ~uplims  # low, high
+            # This is like
+            #     low, high = np.broadcast_to(...)
+            #     return data - low * ~lolims, data + high * ~uplims
+            # except that broadcast_to would strip units.
+            return data + np.row_stack([-(1 - lolims), 1 - uplims]) * err
 
         if xerr is not None:
             left, right = extract_err('x', xerr, x, xlolims, xuplims)

but if you can help with writing a minimal test case to check that, that would be great :-)

Here is a minimal test case from the broken astropy test, though I am not sure how to do it in a way that does not introduce astropy dependency in matplotlib if you are thinking about adding it as a unit test here...

import matplotlib.pyplot as plt

from astropy import units as u
from astropy.visualization.units import quantity_support

with quantity_support():
    x = [1, 2, 3] * u.s
    y = [1, 2, 3] * u.m
    yerr = [3, 2, 1] * u.cm

    fig, ax = plt.subplots()
    ax.errorbar(x, y, yerr=yerr)

    assert ax.xaxis.get_units() == u.s
    assert ax.yaxis.get_units() == u.m

    plt.show()  # Only if you want to see it, not for CI

p.s. Here is the actual implementation of quantity_support -- https://github.com/astropy/astropy/blob/master/astropy/visualization/units.py

@pllim I'm loving that astropy is testing against our default branch :)