Repository: astropy/astropy

Issue:
np.insert gives UnitConversionError with Numpy 1.17.0 (used to work)
Here is a minimal reproducible code snippet:
```python
import numpy as np
from astropy import units as u
PHOTLAM = u.ph / u.cm / u.cm / u.s / u.AA
flux = np.arange(5) * PHOTLAM
newflux = np.insert(flux, 0, 0)
```

With `astropy 4.0.dev25566` and `numpy 1.17.0`, I get this traceback:
```
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
.../astropy/units/quantity.py in to_value(self, unit, equivalencies)
    716             try:
--> 717                 scale = self.unit._to(unit)
    718             except Exception:

.../astropy/units/core.py in _to(self, other)
    952         raise UnitConversionError(
--> 953             f"'{self!r}' is not a scaled version of '{other!r}'")
    954 

UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("ph / (Angstrom cm2 s)")'

During handling of the above exception, another exception occurred:

UnitConversionError                       Traceback (most recent call last)
<ipython-input-5-7e3bcc47515f> in <module>
----> 1 newflux = np.insert(flux, 0, 0)

<__array_function__ internals> in insert(*args, **kwargs)

.../astropy/units/quantity.py in __array_function__(self, function, types, args, kwargs)
   1505             function_helper = FUNCTION_HELPERS[function]
   1506             try:
-> 1507                 args, kwargs, unit, out = function_helper(*args, **kwargs)
   1508             except NotImplementedError:
   1509                 return self._not_implemented_or_raise(function, types)

.../astropy/units/quantity_helper/function_helpers.py in insert(arr, obj, values, *args, **kwargs)
    449         raise NotImplementedError
    450 
--> 451     (arr, values), unit = _quantities2arrays(arr, values)
    452     return (arr, obj, values) + args, kwargs, unit, None
    453 

.../astropy/units/quantity_helper/function_helpers.py in _quantities2arrays(*args)
    319     unit = qs[0].unit
    320     # Allow any units error to be raised.
--> 321     arrays = tuple(q.to_value(unit) for q in qs)
    322     return arrays, unit
    323 

.../astropy/units/quantity_helper/function_helpers.py in <genexpr>(.0)
    319     unit = qs[0].unit
    320     # Allow any units error to be raised.
--> 321     arrays = tuple(q.to_value(unit) for q in qs)
    322     return arrays, unit
    323 

.../astropy/units/quantity.py in to_value(self, unit, equivalencies)
    718             except Exception:
    719                 # Short-cut failed; try default (maybe equivalencies help).
--> 720                 value = self._to_value(unit, equivalencies)
    721             else:
    722                 value = self.view(np.ndarray)

.../astropy/units/quantity.py in _to_value(self, unit, equivalencies)
    649             equivalencies = self._equivalencies
    650         return self.unit.to(unit, self.view(np.ndarray),
--> 651                             equivalencies=equivalencies)
    652 
    653     def to(self, unit, equivalencies=[]):

.../astropy/units/core.py in to(self, other, value, equivalencies)
    987             return UNITY
    988         else:
--> 989             return self._get_converter(other, equivalencies=equivalencies)(value)
    990 
    991     def in_units(self, other, value=1.0, equivalencies=[]):

.../astropy/units/core.py in _get_converter(self, other, equivalencies)
    918                             pass
    919 
--> 920             raise exc
    921 
    922     def _to(self, other):

.../astropy/units/core.py in _get_converter(self, other, equivalencies)
    904         try:
    905             return self._apply_equivalencies(
--> 906                 self, other, self._normalize_equivalencies(equivalencies))
    907         except UnitsError as exc:
    908             # Last hope: maybe other knows how to do it?

.../astropy/units/core.py in _apply_equivalencies(self, unit, other, equivalencies)
    888         raise UnitConversionError(
    889             "{} and {} are not convertible".format(
--> 890                 unit_str, other_str))
    891 
    892     def _get_converter(self, other, equivalencies=[]):

UnitConversionError: '' (dimensionless) and 'ph / (Angstrom cm2 s)' (photon flux density wav) are not convertible
```

With `astropy 4.0.dev25566` and `numpy 1.16.4`, I get this:
```python
>>> newflux  # Another zero inserted at the beginning, as expected
<Quantity [0., 0., 1., 2., 3., 4.] ph / (Angstrom cm2 s)>
```

If I downgrade to `astropy 3.2.1`, it works with both versions of Numpy. I suspect this is a side effect of #8808 . Any advice, @mhvk ?

**UPDATE:** Workaround is to call `np.insert(flux, 0, 0 * flux.unit)` instead.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 5894060123..3bfecd1597 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -15,8 +15,6 @@ astropy.constants
 
 - Default constants now use CODATA 2018 and IAU 2015 definitions. [#8761]
 
-- Constants can be pickled and unpickled. [#9377]
-
 astropy.convolution
 ^^^^^^^^^^^^^^^^^^^
 
diff --git a/astropy/constants/constant.py b/astropy/constants/constant.py
index 3753d1326b..a715cb2fc6 100644
--- a/astropy/constants/constant.py
+++ b/astropy/constants/constant.py
@@ -62,8 +62,7 @@ class ConstantMeta(type):
         exclude = set(['__new__', '__array_finalize__', '__array_wrap__',
                        '__dir__', '__getattr__', '__init__', '__str__',
                        '__repr__', '__hash__', '__iter__', '__getitem__',
-                       '__len__', '__bool__', '__quantity_subclass__',
-                       '__setstate__'])
+                       '__len__', '__bool__', '__quantity_subclass__'])
         for attr, value in vars(Quantity).items():
             if (isinstance(value, types.FunctionType) and
                     attr.startswith('__') and attr.endswith('__') and
diff --git a/astropy/cosmology/core.py b/astropy/cosmology/core.py
index 9af0a8c9d0..c87ee0d98b 100644
--- a/astropy/cosmology/core.py
+++ b/astropy/cosmology/core.py
@@ -1502,8 +1502,10 @@ class FLRW(Cosmology, metaclass=ABCMeta):
           Differential comoving volume per redshift per steradian at
           each input redshift."""
         dh = self._hubble_distance
-        dm = self.comoving_transverse_distance(z)
-        return dh * (dm ** 2.0) / u.Quantity(self.efunc(z), u.steradian)
+        da = self.angular_diameter_distance(z)
+        zp1 = 1.0 + z
+        return dh * ((zp1 * da) ** 2.0) / u.Quantity(self.efunc(z),
+                                                          u.steradian)
 
     def kpc_comoving_per_arcmin(self, z):
         """ Separation in transverse comoving kpc corresponding to an
diff --git a/astropy/units/quantity_helper/function_helpers.py b/astropy/units/quantity_helper/function_helpers.py
index 0a3b5df90a..b708a11cfa 100644
--- a/astropy/units/quantity_helper/function_helpers.py
+++ b/astropy/units/quantity_helper/function_helpers.py
@@ -311,13 +311,42 @@ def _as_quantities(*args):
         raise NotImplementedError
 
 
-def _quantities2arrays(*args):
-    """Convert to Quantities with the unit of the first argument."""
-    qs = _as_quantities(*args)
-    unit = qs[0].unit
-    # Allow any units error to be raised.
-    arrays = tuple(q.to_value(unit) for q in qs)
-    return arrays, unit
+def _quantities2arrays(*args, unit_from_first=False):
+    """Convert to arrays in units of the first argument that has a unit.
+
+    If unit_from_first, take the unit of the first argument regardless
+    whether it actually defined a unit (e.g., dimensionless for arrays).
+    """
+
+    # Turn first argument into a quantity.
+    q = _as_quantity(args[0])
+    if len(args) == 1:
+        return (q.value,), q.unit
+
+    # If we care about the unit being explicit, then check whether this
+    # argument actually had a unit, or was likely inferred.
+    if not unit_from_first and (q.unit is q._default_unit
+                                and not hasattr(args[0], 'unit')):
+        # Here, the argument could still be things like [10*u.one, 11.*u.one]),
+        # i.e., properly dimensionless.  So, we only override with anything
+        # that has a unit not equivalent to dimensionless (fine to ignore other
+        # dimensionless units pass, even if explicitly given).
+        for arg in args[1:]:
+            trial = _as_quantity(arg)
+            if not trial.unit.is_equivalent(q.unit):
+                # Use any explicit unit not equivalent to dimensionless.
+                q = trial
+                break
+
+    # We use the private _to_own_unit method here instead of just
+    # converting everything to quantity and then do .to_value(qs0.unit)
+    # as we want to allow arbitrary unit for 0, inf, and nan.
+    try:
+        arrays = tuple((q._to_own_unit(arg)) for arg in args)
+    except TypeError:
+        raise NotImplementedError
+
+    return arrays, q.unit
 
 
 def _iterable_helper(*args, out=None, **kwargs):
@@ -435,7 +464,7 @@ def piecewise(x, condlist, funclist, *args, **kw):
 
 @function_helper
 def append(arr, values, *args, **kwargs):
-    arrays, unit = _quantities2arrays(arr, values)
+    arrays, unit = _quantities2arrays(arr, values, unit_from_first=True)
     return arrays + args, kwargs, unit, None
 
 
@@ -446,7 +475,8 @@ def insert(arr, obj, values, *args, **kwargs):
     if isinstance(obj, Quantity):
         raise NotImplementedError
 
-    (arr, values), unit = _quantities2arrays(arr, values)
+    (arr, values), unit = _quantities2arrays(arr, values,
+                                             unit_from_first=True)
     return (arr, obj, values) + args, kwargs, unit, None
 
 
@@ -573,7 +603,7 @@ def bincount(x, weights=None, minlength=0):
 
 @function_helper
 def digitize(x, bins, *args, **kwargs):
-    arrays, unit = _quantities2arrays(x, bins)
+    arrays, unit = _quantities2arrays(x, bins, unit_from_first=True)
     return arrays + args, kwargs, None, None
 
 
diff --git a/docs/conf.py b/docs/conf.py
index 386f3bf2fc..bfaf8fd40f 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -192,7 +192,6 @@ man_pages = [('index', project.lower(), project + u' Documentation',
 
 # Setting this URL is requited by sphinx-astropy
 github_issues_url = 'https://github.com/astropy/astropy/issues/'
-edit_on_github_branch = 'master'
 
 # Enable nitpicky mode - which ensures that all references in the docs
 # resolve.


Test patch:
diff --git a/astropy/constants/tests/test_pickle.py b/astropy/constants/tests/test_pickle.py
index 75c2fb93bf..b9714f371f 100644
--- a/astropy/constants/tests/test_pickle.py
+++ b/astropy/constants/tests/test_pickle.py
@@ -9,8 +9,8 @@ originals = [const.Constant('h_fake', 'Not Planck',
                             0.0, 'J s', 0.0, 'fakeref',
                             system='si'),
              const.h,
-             const.e.si]
-xfails = [False, False, False]
+             const.e]
+xfails = [True, True, True]
 
 
 @pytest.mark.parametrize(("original", "xfail"), zip(originals, xfails))
diff --git a/astropy/units/tests/test_quantity_non_ufuncs.py b/astropy/units/tests/test_quantity_non_ufuncs.py
index c71a42d616..af422b241f 100644
--- a/astropy/units/tests/test_quantity_non_ufuncs.py
+++ b/astropy/units/tests/test_quantity_non_ufuncs.py
@@ -436,10 +436,10 @@ class TestConcatenate(metaclass=CoverageMeta):
 
     def check(self, func, *args, **kwargs):
         q_list = kwargs.pop('q_list', [self.q1, self.q2])
+        q_ref = kwargs.pop('q_ref', q_list[0])
         o = func(q_list, *args, **kwargs)
-        unit = q_list[0].unit
-        v_list = [q.to_value(unit) for q in q_list]
-        expected = func(v_list, *args, **kwargs) * unit
+        v_list = [q_ref._to_own_unit(q) for q in q_list]
+        expected = func(v_list, *args, **kwargs) * q_ref.unit
         assert o.shape == expected.shape
         assert np.all(o == expected)
 
@@ -449,6 +449,9 @@ class TestConcatenate(metaclass=CoverageMeta):
         self.check(np.concatenate)
         self.check(np.concatenate, axis=1)
 
+        self.check(np.concatenate, q_list=[np.zeros(self.q1.shape), self.q1, self.q2],
+                   q_ref=self.q1)
+
         out = np.empty((4, 3)) * u.dimensionless_unscaled
         result = np.concatenate([self.q1, self.q2], out=out)
         assert out is result
@@ -490,6 +493,9 @@ class TestConcatenate(metaclass=CoverageMeta):
     def test_block(self):
         self.check(np.block)
 
+        result = np.block([[0., 1.*u.m], [1.*u.cm, 2.*u.km]])
+        assert np.all(result == np.block([[0, 1.], [.01, 2000.]]) << u.m)
+
     @pytest.mark.xfail(NO_ARRAY_FUNCTION,
                        reason="Needs __array_function__ support")
     def test_append(self):
@@ -509,13 +515,17 @@ class TestConcatenate(metaclass=CoverageMeta):
     @pytest.mark.xfail(NO_ARRAY_FUNCTION,
                        reason="Needs __array_function__ support")
     def test_insert(self):
-        # Unit of inserted values is ignored.
+        # Unit of inserted values is not ignored.
         q = np.arange(12.).reshape(6, 2) * u.m
         out = np.insert(q, (3, 5), [50., 25.] * u.cm)
         assert isinstance(out, u.Quantity)
         assert out.unit == q.unit
-        expected = np.insert(q.value, (3, 5), [0.5, 0.25]) * u.m
+        expected = np.insert(q.value, (3, 5), [0.5, 0.25]) << q.unit
         assert np.all(out == expected)
+        # 0 can have any unit.
+        out2 = np.insert(q, (3, 5), 0)
+        expected2 = np.insert(q.value, (3, 5), 0) << q.unit
+        assert np.all(out2 == expected2)
 
         a = np.arange(3.)
         result = np.insert(a, (2,), 50. * u.percent)
@@ -526,6 +536,8 @@ class TestConcatenate(metaclass=CoverageMeta):
 
         with pytest.raises(TypeError):
             np.insert(q, 3 * u.cm, 50. * u.cm)
+        with pytest.raises(u.UnitsError):
+            np.insert(q, (3, 5), 0. * u.s)
 
     @pytest.mark.xfail(NO_ARRAY_FUNCTION,
                        reason="Needs __array_function__ support")
@@ -1211,6 +1223,13 @@ class TestInterpolationFunctions(metaclass=CoverageMeta):
         assert out2.unit == expected2.unit
         assert np.all(out2 == expected2)
 
+        out3 = np.piecewise(x, [x < 1 * u.m, x >= 0],
+                            [0, 1*u.percent, lambda x: 1*u.one])
+        expected3 = np.piecewise(x.value, [x.value < 1, x.value >= 0],
+                                 [0, 0.01, 1]) * u.one
+        assert out3.unit == expected3.unit
+        assert np.all(out3 == expected3)
+
         with pytest.raises(TypeError):  # no Quantity in condlist.
             np.piecewise(x, [x], [0.])
 
@@ -1692,6 +1711,9 @@ class TestSetOpsFcuntions(metaclass=CoverageMeta):
                        reason="Needs __array_function__ support")
     def test_union1d(self):
         self.check2(np.union1d)
+        result = np.union1d(np.array([0., np.nan]), np.arange(3) << u.m)
+        assert result.unit is u.m
+        assert_array_equal(result.value, np.array([0., 1., 2., np.nan]))
 
     @pytest.mark.xfail(NO_ARRAY_FUNCTION,
                        reason="Needs __array_function__ support")
@@ -1702,6 +1724,10 @@ class TestSetOpsFcuntions(metaclass=CoverageMeta):
                        reason="Needs __array_function__ support")
     def test_in1d(self):
         self.check2(np.in1d, unit=None)
+        # Check zero is treated as having any unit.
+        assert np.in1d(np.zeros(1), self.q2)
+        with pytest.raises(u.UnitsError):
+            np.in1d(np.ones(1), self.q2)
 
     @pytest.mark.xfail(NO_ARRAY_FUNCTION,
                        reason="Needs __array_function__ support")