np.median drops units if the result is nan

[mhvk]

As reported by @larrybradley in #12146 (comment):
np.median drops the unit on Quantity input where the result is nan:

>>> import numpy as np
>>> import astropy.units as u
>>> data = np.array([1., 2, np.nan, 3, 4]) << u.km
>>> np.median(data)
nan
>>> np.mean(data)
<Quantity nan km>
>>> np.var(data)
<Quantity nan km2>

Note that it works fine on multi-dimensional output:

In [6]: np.median([[1, 2], [np.nan, 3]] << u.km, axis=1)
Out[6]: <Quantity [1.5, nan] km>
In [8]: np.median([[1, np.nan], [np.nan, 3]] << u.km, axis=1)
Out[8]: <Quantity [nan, nan] km>

[mhvk]

At some level, this reflects a problem upstream: while generally the numpy code deals well with subclasses, it seems it does not for this particular case - which I clearly had not noticed, so np.median is listed as a "subclass-safe" function at

astropy/astropy/units/quantity_helper/function_helpers.py

Line 82 in 23404ee

np.average, np.mean, np.std, np.var, np.median, np.trace,

[nstarman]

@mhvk, is this something we could fix with #11893 ?

I'm planning on returning to this soon as my PR for to/from_format("astropy.model") requires a quantity-aware np.vectorized.

[mhvk]

@nstarman - no, it is really a bit of a numpy bug, which we can work around in function_helpers.

[mhvk]

See numpy/numpy#19869 (review) for the upstream fix. We could fix it for existing numpy versions if need be.

[nstarman]

Would that be in compat or a median-specific ufunc override?

[mhvk]

Would that be in compat or a median-specific ufunc override?

I was thinking of a numpy-version dependent override of np.median in units.quantity_helpers.function_helpers. Obviously, if we go that route, the work-arounds in stats can be removed.