Repository: astropy/astropy

Issue:
Primary HDU checksum data must be np.array?
I have some code which attempts to update the FITS checksums:

     output = "out.fits"
     with fits.open(self.filepath, do_not_scale_image_data=True) as hdus:
             for hdu in hdus:
                   fits.append(output, hdu.data, hdu.header, checksum=True)

For a file which has no data in an HDU (primary),  this code now raises an exception requiring that the data array be an np.array rather than None.

There's a relatively easy "work" around something like this:

         data = hdu.data if hdu.data is not None else np.array([])

but I wasn't positive it's correct.

Since it used to be possible to round trip these HDUs,  I thought I should raise it as an issue and possible bug.

I (think) this problem is exposed with:

astropy                   4.0.0.dev12788+g4166b6a94          py36_0    http://ssb.stsci.edu/conda-dev
numpy                     1.16.3           py36hacdab7b_0  
numpy-base                1.16.3           py36h6575580_0  



Discussion:
pllim: cc @mcara 
MSeifert04: Puh, these are those kind of consequences that make `io.fits` so "interesting" 😕 

I don't have the time atm to do a proper in-depth review. However if it fixes the reported problem and the test suite doesn't complain it seems like a clear improvement.
saimn: I could also revert the change to `_BaseHDU.__new__`, but given that no test was failing because of this, I think it's a good point to avoid the double detection from the header, and introduce this classmethod for specific cases like this one. Hopefully not having the detection in `_BaseHDU.__new__` should not have other unintended consequences :crossed_fingers: 
pllim: @jaytmiller , can you please confirm that works for you?
pllim: Well, that's weird. There is a `html-docs` error that I don't see on latest `master`. Maybe rebasing would help? 

```
/home/circleci/project/docs/whatsnew/3.2.rst:98:Inline literal start-string without end-string.
```
bsipocz: I believe that failure got fixed by #8669
astrofrog: @saimn - can you rebase? Thanks!
codecov[bot]: # [Codecov](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=h1) Report
> Merging [#8666](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=desc) into [master](https://codecov.io/gh/astropy/astropy/commit/8a59b161fc83957bbc39efafcfbf4ad36b973373?src=pr&el=desc) will **decrease** coverage by `<.01%`.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/astropy/astropy/pull/8666/graphs/tree.svg?width=650&token=eEhbJxwU0j&height=150&src=pr)](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #8666      +/-   ##
==========================================
- Coverage   86.95%   86.95%   -0.01%     
==========================================
  Files         399      399              
  Lines       59351    59335      -16     
  Branches     1100     1100              
==========================================
- Hits        51609    51594      -15     
+ Misses       7101     7100       -1     
  Partials      641      641
```


| [Impacted Files](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [astropy/io/fits/convenience.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS9pby9maXRzL2NvbnZlbmllbmNlLnB5) | `85.76% <100%> (ø)` | :arrow_up: |
| [astropy/io/fits/hdu/base.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS9pby9maXRzL2hkdS9iYXNlLnB5) | `90.36% <100%> (+0.04%)` | :arrow_up: |
| [astropy/wcs/wcsapi/fitswcs.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS93Y3Mvd2NzYXBpL2ZpdHN3Y3MucHk=) | `91.91% <0%> (-0.18%)` | :arrow_down: |
| [astropy/wcs/wcsapi/high\_level\_api.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS93Y3Mvd2NzYXBpL2hpZ2hfbGV2ZWxfYXBpLnB5) | `97.97% <0%> (-0.08%)` | :arrow_down: |
| [astropy/io/fits/header.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS9pby9maXRzL2hlYWRlci5weQ==) | `96% <0%> (-0.05%)` | :arrow_down: |
| [astropy/io/fits/card.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS9pby9maXRzL2NhcmQucHk=) | `86.47% <0%> (-0.05%)` | :arrow_down: |
| [astropy/wcs/wcsapi/low\_level\_api.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS93Y3Mvd2NzYXBpL2xvd19sZXZlbF9hcGkucHk=) | `97.05% <0%> (ø)` | :arrow_up: |
| [astropy/io/fits/hdu/image.py](https://codecov.io/gh/astropy/astropy/pull/8666/diff?src=pr&el=tree#diff-YXN0cm9weS9pby9maXRzL2hkdS9pbWFnZS5weQ==) | `96.25% <0%> (+0.2%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=footer). Last update [8a59b16...d665d91](https://codecov.io/gh/astropy/astropy/pull/8666?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

saimn: @astrofrog - Done.
saimn: Anyone to review/approve? :pray: 
jaytmiller: Sorry,  I looked this over and was relatively convinced it should work based on the  regression test case.
pllim: Doc build is still failing after rebase and its log indicates that it went over the max 10 minutes limit.

```
...
reading sources... [ 95%] timeseries/times
reading sources... [ 95%] uncertainty/index
Too long with no output (exceeded 10m0s)
```

Does this PR somehow make the doc build too long? @astrofrog or @Cadair , how do we increase the limit on CircleCI?
saimn: This PR should not increase the build time, and the error looks similar to the ones with the recent mpl issue. I can try a new rebase to be sure?
pllim: > I can try a new rebase to be sure?

Feel free to try again. Thanks!

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 1026ed63b1..13bf5e6d75 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -311,10 +311,6 @@ astropy.uncertainty
 astropy.units
 ^^^^^^^^^^^^^
 
-- Add support for the ``clip`` ufunc, which in numpy 1.17 is used to implement
-  ``np.clip``.  As part of that, remove the ``Quantity.clip`` method under
-  numpy 1.17. [#8747]
-
 astropy.utils
 ^^^^^^^^^^^^^
 
diff --git a/astropy/io/fits/convenience.py b/astropy/io/fits/convenience.py
index 84d30d2b64..ac5c25c5ff 100644
--- a/astropy/io/fits/convenience.py
+++ b/astropy/io/fits/convenience.py
@@ -1034,7 +1034,7 @@ def _getext(filename, mode, *args, ext=None, extname=None, extver=None,
 def _makehdu(data, header):
     if header is None:
         header = Header()
-    hdu = _BaseHDU(data, header)
+    hdu = _BaseHDU._from_data(data, header)
     if hdu.__class__ in (_BaseHDU, _ValidHDU):
         # The HDU type was unrecognized, possibly due to a
         # nonexistent/incomplete header
diff --git a/astropy/io/fits/hdu/base.py b/astropy/io/fits/hdu/base.py
index accc72ee30..6f8a0b2457 100644
--- a/astropy/io/fits/hdu/base.py
+++ b/astropy/io/fits/hdu/base.py
@@ -56,8 +56,19 @@ class InvalidHDUException(Exception):
 
 def _hdu_class_from_header(cls, header):
     """
-    Used primarily by _BaseHDU.__new__ to find an appropriate HDU class to use
-    based on values in the header.  See the _BaseHDU.__new__ docstring.
+    Iterates through the subclasses of _BaseHDU and uses that class's
+    match_header() method to determine which subclass to instantiate.
+
+    It's important to be aware that the class hierarchy is traversed in a
+    depth-last order.  Each match_header() should identify an HDU type as
+    uniquely as possible.  Abstract types may choose to simply return False
+    or raise NotImplementedError to be skipped.
+
+    If any unexpected exceptions are raised while evaluating
+    match_header(), the type is taken to be _CorruptedHDU.
+
+    Used primarily by _BaseHDU._readfrom_internal and _BaseHDU._from_data to
+    find an appropriate HDU class to use based on values in the header.
     """
 
     klass = cls  # By default, if no subclasses are defined
@@ -361,6 +372,15 @@ class _BaseHDU(metaclass=_BaseHDUMeta):
         hdulist.writeto(name, output_verify, overwrite=overwrite,
                         checksum=checksum)
 
+    @classmethod
+    def _from_data(cls, data, header, **kwargs):
+        """
+        Instantiate the HDU object after guessing the HDU class from the
+        FITS Header.
+        """
+        klass = _hdu_class_from_header(cls, header)
+        return klass(data=data, header=header, **kwargs)
+
     @classmethod
     def _readfrom_internal(cls, data, header=None, checksum=False,
                            ignore_missing_end=False, **kwargs):
diff --git a/astropy/units/function/core.py b/astropy/units/function/core.py
index 50cca5bf64..0195770474 100644
--- a/astropy/units/function/core.py
+++ b/astropy/units/function/core.py
@@ -683,7 +683,3 @@ class FunctionQuantity(Quantity):
 
     def cumsum(self, axis=None, dtype=None, out=None):
         return self._wrap_function(np.cumsum, axis, dtype, out=out)
-
-    def clip(self, a_min, a_max, out=None):
-        return self._wrap_function(np.clip, self._to_own_unit(a_min),
-                                   self._to_own_unit(a_max), out=out)
diff --git a/astropy/units/quantity.py b/astropy/units/quantity.py
index da5d7e0a55..73c5c66e69 100644
--- a/astropy/units/quantity.py
+++ b/astropy/units/quantity.py
@@ -20,7 +20,7 @@ from .core import (Unit, dimensionless_unscaled, get_current_unit_registry,
                    UnitBase, UnitsError, UnitConversionError, UnitTypeError)
 from .utils import is_effectively_unity
 from .format.latex import Latex
-from astropy.utils.compat import NUMPY_LT_1_14, NUMPY_LT_1_16, NUMPY_LT_1_17
+from astropy.utils.compat import NUMPY_LT_1_14, NUMPY_LT_1_16
 from astropy.utils.compat.misc import override__dir__
 from astropy.utils.exceptions import AstropyDeprecationWarning, AstropyWarning
 from astropy.utils.misc import isiterable, InheritDocstrings
@@ -455,10 +455,9 @@ class Quantity(np.ndarray, metaclass=InheritDocstrings):
             kwargs['out'] = (out_array,) if function.nout == 1 else out_array
 
         # Same for inputs, but here also convert if necessary.
-        arrays = []
-        for input_, converter in zip(inputs, converters):
-            input_ = getattr(input_, 'value', input_)
-            arrays.append(converter(input_) if converter else input_)
+        arrays = [(converter(input_.value) if converter else
+                   getattr(input_, 'value', input_))
+                  for input_, converter in zip(inputs, converters)]
 
         # Call our superclass's __array_ufunc__
         result = super().__array_ufunc__(function, method, *arrays, **kwargs)
@@ -1502,10 +1501,9 @@ class Quantity(np.ndarray, metaclass=InheritDocstrings):
         result = function(*args, **kwargs)
         return self._result_as_quantity(result, unit, out)
 
-    if NUMPY_LT_1_17:
-        def clip(self, a_min, a_max, out=None):
-            return self._wrap_function(np.clip, self._to_own_unit(a_min),
-                                       self._to_own_unit(a_max), out=out)
+    def clip(self, a_min, a_max, out=None):
+        return self._wrap_function(np.clip, self._to_own_unit(a_min),
+                                   self._to_own_unit(a_max), out=out)
 
     def trace(self, offset=0, axis1=0, axis2=1, dtype=None, out=None):
         return self._wrap_function(np.trace, offset, axis1, axis2, dtype,
diff --git a/astropy/units/quantity_helper/converters.py b/astropy/units/quantity_helper/converters.py
index cfcd744589..570655bfcb 100644
--- a/astropy/units/quantity_helper/converters.py
+++ b/astropy/units/quantity_helper/converters.py
@@ -166,34 +166,31 @@ def converters_and_unit(function, method, *args):
         converters, result_unit = ufunc_helper(function, *units)
 
         if any(converter is False for converter in converters):
-            # for multi-argument ufuncs with a quantity and a non-quantity,
+            # for two-argument ufuncs with a quantity and a non-quantity,
             # the quantity normally needs to be dimensionless, *except*
             # if the non-quantity can have arbitrary unit, i.e., when it
             # is all zero, infinity or NaN.  In that case, the non-quantity
             # can just have the unit of the quantity
             # (this allows, e.g., `q > 0.` independent of unit)
+            maybe_arbitrary_arg = args[converters.index(False)]
             try:
-                # Don't fold this loop in the test above: this rare case
-                # should not make the common case slower.
-                for i, converter in enumerate(converters):
-                    if converter is not False:
-                        continue
-                    if can_have_arbitrary_unit(args[i]):
-                        converters[i] = None
-                    else:
-                        raise UnitConversionError(
-                            "Can only apply '{0}' function to "
-                            "dimensionless quantities when other "
-                            "argument is not a quantity (unless the "
-                            "latter is all zero/infinity/nan)"
-                            .format(function.__name__))
+                if can_have_arbitrary_unit(maybe_arbitrary_arg):
+                    converters = [None, None]
+                else:
+                    raise UnitConversionError(
+                        "Can only apply '{0}' function to "
+                        "dimensionless quantities when other "
+                        "argument is not a quantity (unless the "
+                        "latter is all zero/infinity/nan)"
+                        .format(function.__name__))
             except TypeError:
                 # _can_have_arbitrary_unit failed: arg could not be compared
                 # with zero or checked to be finite. Then, ufunc will fail too.
                 raise TypeError("Unsupported operand type(s) for ufunc {0}: "
-                                "'{1}'".format(function.__name__,
-                                               ','.join([arg.__class__.__name__
-                                                         for arg in args])))
+                                "'{1}' and '{2}'"
+                                .format(function.__name__,
+                                        args[0].__class__.__name__,
+                                        args[1].__class__.__name__))
 
         # In the case of np.power and np.float_power, the unit itself needs to
         # be modified by an amount that depends on one of the input values,
diff --git a/astropy/units/quantity_helper/helpers.py b/astropy/units/quantity_helper/helpers.py
index d08c69352f..5bbd4d0e68 100644
--- a/astropy/units/quantity_helper/helpers.py
+++ b/astropy/units/quantity_helper/helpers.py
@@ -296,39 +296,6 @@ def helper_divmod(f, unit1, unit2):
     return converters, (dimensionless_unscaled, result_unit)
 
 
-def helper_clip(f, unit1, unit2, unit3):
-    # Treat the array being clipped as primary.
-    converters = [None]
-    if unit1 is None:
-        result_unit = dimensionless_unscaled
-        try:
-            converters += [(None if unit is None else
-                            get_converter(unit, dimensionless_unscaled))
-                           for unit in (unit2, unit3)]
-        except UnitsError:
-            raise UnitConversionError(
-                "Can only apply '{0}' function to quantities with "
-                "compatible dimensions".format(f.__name__))
-
-    else:
-        result_unit = unit1
-        for unit in unit2, unit3:
-            try:
-                converter = get_converter(_d(unit), result_unit)
-            except UnitsError:
-                if unit is None:
-                    # special case: OK if unitless number is zero, inf, nan
-                    converters.append(False)
-                else:
-                    raise UnitConversionError(
-                        "Can only apply '{0}' function to quantities with "
-                        "compatible dimensions".format(f.__name__))
-            else:
-                converters.append(converter)
-
-    return converters, result_unit
-
-
 # list of ufuncs:
 # http://docs.scipy.org/doc/numpy/reference/ufuncs.html#available-ufuncs
 
@@ -440,6 +407,3 @@ UFUNC_HELPERS[np.floor_divide] = helper_twoarg_floor_divide
 UFUNC_HELPERS[np.heaviside] = helper_heaviside
 UFUNC_HELPERS[np.float_power] = helper_power
 UFUNC_HELPERS[np.divmod] = helper_divmod
-# Check for clip ufunc; note that np.clip is a wrapper function, not the ufunc.
-if isinstance(getattr(np.core.umath, 'clip', None), np.ufunc):
-    UFUNC_HELPERS[np.core.umath.clip] = helper_clip
diff --git a/astropy/utils/compat/numpycompat.py b/astropy/utils/compat/numpycompat.py
index 7db074894f..23d3a4d4ac 100644
--- a/astropy/utils/compat/numpycompat.py
+++ b/astropy/utils/compat/numpycompat.py
@@ -7,7 +7,7 @@ from astropy.utils import minversion
 
 
 __all__ = ['NUMPY_LT_1_14', 'NUMPY_LT_1_14_1', 'NUMPY_LT_1_14_2',
-           'NUMPY_LT_1_16', 'NUMPY_LT_1_17']
+           'NUMPY_LT_1_16']
 
 # TODO: It might also be nice to have aliases to these named for specific
 # features/bugs we're checking for (ex:
@@ -16,4 +16,3 @@ NUMPY_LT_1_14 = not minversion('numpy', '1.14')
 NUMPY_LT_1_14_1 = not minversion('numpy', '1.14.1')
 NUMPY_LT_1_14_2 = not minversion('numpy', '1.14.2')
 NUMPY_LT_1_16 = not minversion('numpy', '1.16')
-NUMPY_LT_1_17 = not minversion('numpy', '1.17')


Test patch:
diff --git a/astropy/io/fits/tests/test_convenience.py b/astropy/io/fits/tests/test_convenience.py
index f0702d1742..7f475291e7 100644
--- a/astropy/io/fits/tests/test_convenience.py
+++ b/astropy/io/fits/tests/test_convenience.py
@@ -187,3 +187,17 @@ class TestConvenience(FitsTestCase):
         append_file = tmpdir.join('append.fits')
         with append_file.open(mode) as handle:
             fits.append(filename=handle, data=np.ones((4, 4)))
+
+    def test_append_with_header(self):
+        """
+        Test fits.append with a fits Header, which triggers detection of the
+        HDU class. Regression test for
+        https://github.com/astropy/astropy/issues/8660
+        """
+        testfile = self.temp('test_append_1.fits')
+        with fits.open(self.data('test0.fits')) as hdus:
+            for hdu in hdus:
+                fits.append(testfile, hdu.data, hdu.header, checksum=True)
+
+        with fits.open(testfile, checksum=True) as hdus:
+            assert len(hdus) == 5
diff --git a/astropy/units/tests/test_quantity.py b/astropy/units/tests/test_quantity.py
index 9f353a096d..13cb6b4d1f 100644
--- a/astropy/units/tests/test_quantity.py
+++ b/astropy/units/tests/test_quantity.py
@@ -496,10 +496,11 @@ class TestQuantityOperations:
 
     def test_non_number_type(self):
         q1 = u.Quantity(11.412, unit=u.meter)
+        type_err_msg = ("Unsupported operand type(s) for ufunc add: "
+                        "'Quantity' and 'dict'")
         with pytest.raises(TypeError) as exc:
             q1 + {'a': 1}
-        assert exc.value.args[0].startswith(
-            "Unsupported operand type(s) for ufunc add:")
+        assert exc.value.args[0] == type_err_msg
 
         with pytest.raises(TypeError):
             q1 + u.meter
diff --git a/astropy/units/tests/test_quantity_ufuncs.py b/astropy/units/tests/test_quantity_ufuncs.py
index ac0db53374..49ffe4e907 100644
--- a/astropy/units/tests/test_quantity_ufuncs.py
+++ b/astropy/units/tests/test_quantity_ufuncs.py
@@ -868,92 +868,6 @@ class TestInplaceUfuncs:
             a4 += u.Quantity(10, u.mm, dtype=np.int64)
 
 
-@pytest.mark.skipif(not hasattr(np.core.umath, 'clip'),
-                    reason='no clip ufunc available')
-class TestClip:
-    """Test the clip ufunc.
-
-    In numpy, this is hidden behind a function that does not backwards
-    compatibility checks.  We explicitly test the ufunc here.
-    """
-    def setup(self):
-        self.clip = np.core.umath.clip
-
-    def test_clip_simple(self):
-        q = np.arange(-1., 10.) * u.m
-        q_min = 125 * u.cm
-        q_max = 0.0055 * u.km
-        result = self.clip(q, q_min, q_max)
-        assert result.unit == q.unit
-        expected = self.clip(q.value, q_min.to_value(q.unit),
-                             q_max.to_value(q.unit)) * q.unit
-        assert np.all(result == expected)
-
-    def test_clip_unitless_parts(self):
-        q = np.arange(-1., 10.) * u.m
-        qlim = 0.0055 * u.km
-        # one-sided
-        result1 = self.clip(q, -np.inf, qlim)
-        expected1 = self.clip(q.value, -np.inf, qlim.to_value(q.unit)) * q.unit
-        assert np.all(result1 == expected1)
-        result2 = self.clip(q, qlim, np.inf)
-        expected2 = self.clip(q.value, qlim.to_value(q.unit), np.inf) * q.unit
-        assert np.all(result2 == expected2)
-        # Zero
-        result3 = self.clip(q, np.zeros(q.shape), qlim)
-        expected3 = self.clip(q.value, 0, qlim.to_value(q.unit)) * q.unit
-        assert np.all(result3 == expected3)
-        # Two unitless parts, array-shaped.
-        result4 = self.clip(q, np.zeros(q.shape), np.full(q.shape, np.inf))
-        expected4 = self.clip(q.value, 0, np.inf) * q.unit
-        assert np.all(result4 == expected4)
-
-    def test_clip_dimensionless(self):
-        q = np.arange(-1., 10.) * u.dimensionless_unscaled
-        result = self.clip(q, 200 * u.percent, 5.)
-        expected = self.clip(q, 2., 5.)
-        assert result.unit == u.dimensionless_unscaled
-        assert np.all(result == expected)
-
-    def test_clip_ndarray(self):
-        a = np.arange(-1., 10.)
-        result = self.clip(a, 200 * u.percent, 5. * u.dimensionless_unscaled)
-        assert isinstance(result, u.Quantity)
-        expected = self.clip(a, 2., 5.) * u.dimensionless_unscaled
-        assert np.all(result == expected)
-
-    def test_clip_quantity_inplace(self):
-        q = np.arange(-1., 10.) * u.m
-        q_min = 125 * u.cm
-        q_max = 0.0055 * u.km
-        expected = self.clip(q.value, q_min.to_value(q.unit),
-                             q_max.to_value(q.unit)) * q.unit
-        result = self.clip(q, q_min, q_max, out=q)
-        assert result is q
-        assert np.all(result == expected)
-
-    def test_clip_ndarray_dimensionless_output(self):
-        a = np.arange(-1., 10.)
-        q = np.zeros_like(a) * u.m
-        expected = self.clip(a, 2., 5.) * u.dimensionless_unscaled
-        result = self.clip(a, 200 * u.percent, 5. * u.dimensionless_unscaled,
-                           out=q)
-        assert result is q
-        assert result.unit == u.dimensionless_unscaled
-        assert np.all(result == expected)
-
-    def test_clip_errors(self):
-        q = np.arange(-1., 10.) * u.m
-        with pytest.raises(u.UnitsError):
-            self.clip(q, 0, 1*u.s)
-        with pytest.raises(u.UnitsError):
-            self.clip(q.value, 0, 1*u.s)
-        with pytest.raises(u.UnitsError):
-            self.clip(q, -1, 0.)
-        with pytest.raises(u.UnitsError):
-            self.clip(q, 0., 1.)
-
-
 class TestUfuncAt:
     """Test that 'at' method for ufuncs (calculates in-place at given indices)