wcs.world_axis_physical_types check breaks apStar_loader

[dhomeier]

Since merging #754 its tests fail as

________________________________________________________________________________ test_apStar_loader ________________________________________________________________________________

tmpdir = local('/private/var/folders/bb/f_z2vsjd21d_p0jpqxzhhwhh00027w/T/pytest-of-derek/pytest-52/test_apStar_loader0')

    @pytest.mark.remote_data
    def test_apStar_loader(tmpdir):
        apstar_url = ("https://data.sdss.org/sas/dr16/apogee/spectro/redux/r12/"
                      "stars/apo25m/N7789/apStar-r12-2M00005414+5522241.fits")
        with set_temp_cache(path=str(tmpdir)):
            filename = download_file(apstar_url, cache=True)
>           spectrum = apStar_loader(filename)

lib/python3.8/site-packages/specutils/io/default_loaders/tests/test_apogee.py:17: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
lib/python3.8/site-packages/specutils/io/registers.py:111: in wrapper
    return func(*args, **kwargs)
lib/python3.8/site-packages/specutils/io/default_loaders/apogee.py:146: in apStar_loader
    return Spectrum1D(data=data * unit,
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <[AttributeError("'Spectrum1D' object has no attribute '_data'") raised in repr()] Spectrum1D object at 0x12cb64fd0>
flux = <Quantity [0., 0., 0., ..., 0., 0., 0.] 1e-17 erg / (Angstrom cm2 s)>
spectral_axis = <Quantity [15100.80154164, 15101.01016837, 15101.21879797, ...,
           16999.33764338, 16999.57249953, 16999.80735892] Angstrom>
wcs = WCS Keywords

Number of WCS axes: 2
CTYPE : 'LOG-LINEAR'  ''  
CRVAL : 4.179  0.0  
CRPIX : 1.0  0.0  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : 6e-06  1.0  
NAXIS : 8575  9
velocity_convention = None, rest_value = None, redshift = None, radial_velocity = None, bin_specification = 'centers'
kwargs = {'meta': {'header': SIMPLE  =                    T /image conforms to FITS standard                 
BITPIX  =        ...StdDevUncertainty([1.e+10, 1.e+10, 1.e+10, ..., 1.e+10, 1.e+10, 1.e+10]), 'unit': Unit("1e-17 erg / (Angstrom cm2 s)")}
unknown_kwargs = set(), temp_axes = [], phys_axes = [None, None], i = 0

    def __init__(self, flux=None, spectral_axis=None, wcs=None,
                 velocity_convention=None, rest_value=None, redshift=None,
                 radial_velocity=None, bin_specification=None, **kwargs):
        # Check for pre-defined entries in the kwargs dictionary.
        unknown_kwargs = set(kwargs).difference(
            {'data', 'unit', 'uncertainty', 'meta', 'mask', 'copy'})
    
        if len(unknown_kwargs) > 0:
            raise ValueError("Initializer contains unknown arguments(s): {}."
                             "".format(', '.join(map(str, unknown_kwargs))))
    
        # If the flux (data) argument is a subclass of nddataref (as it would
        # be for internal arithmetic operations), avoid setup entirely.
        if isinstance(flux, NDDataRef):
            super().__init__(flux)
            return
    
        # If the mask kwarg is not passed to the constructor, but the flux array
        # contains NaNs, add the NaN locations to the mask.
        if "mask" not in kwargs and flux is not None:
            nan_mask = np.isnan(flux)
            if nan_mask.any():
                if hasattr(self, "mask"):
                    kwargs["mask"] = np.logical_or(nan_mask, self.mask)
                else:
                    kwargs["mask"] = nan_mask.copy()
            del nan_mask
    
        # Ensure that the flux argument is an astropy quantity
        if flux is not None:
            if not isinstance(flux, u.Quantity):
                raise ValueError("Flux must be a `Quantity` object.")
            elif flux.isscalar:
                flux = u.Quantity([flux])
    
        # Ensure that only one or neither of these parameters is set
        if redshift is not None and radial_velocity is not None:
            raise ValueError("Cannot set both radial_velocity and redshift at "
                             "the same time.")
    
        # In cases of slicing, new objects will be initialized with `data`
        # instead of ``flux``. Ensure we grab the `data` argument.
        if flux is None and 'data' in kwargs:
            flux = kwargs.pop('data')
    
        # Ensure that the unit information codified in the quantity object is
        # the One True Unit.
        kwargs.setdefault('unit', flux.unit if isinstance(flux, u.Quantity)
                                            else kwargs.get('unit'))
    
        # In the case where the arithmetic operation is being performed with
        # a single float, int, or array object, just go ahead and ignore wcs
        # requirements
        if (not isinstance(flux, u.Quantity) or isinstance(flux, float)
            or isinstance(flux, int)) and np.ndim(flux) == 0:
    
            super(Spectrum1D, self).__init__(data=flux, wcs=wcs, **kwargs)
            return
    
        if rest_value is None:
            if hasattr(wcs, 'rest_frequency') and wcs.rest_frequency != 0:
                rest_value = wcs.rest_frequency * u.Hz
            elif hasattr(wcs, 'rest_wavelength') and wcs.rest_wavelength != 0:
                rest_value = wcs.rest_wavelength * u.AA
            elif hasattr(wcs, 'wcs') and hasattr(wcs.wcs, 'restfrq') and wcs.wcs.restfrq > 0:
                rest_value = wcs.wcs.restfrq * u.Hz
            elif hasattr(wcs, 'wcs') and hasattr(wcs.wcs, 'restwav') and wcs.wcs.restwav > 0:
                rest_value = wcs.wcs.restwav * u.m
            else:
                rest_value = None
        else:
            if not isinstance(rest_value, u.Quantity):
                log.info("No unit information provided with rest value. "
                             "Assuming units of spectral axis ('%s').",
                             spectral_axis.unit)
                rest_value = u.Quantity(rest_value, spectral_axis.unit)
            elif not rest_value.unit.is_equivalent(u.AA, equivalencies=u.spectral()):
                raise u.UnitsError("Rest value must be "
                                   "energy/wavelength/frequency equivalent.")
    
        # If flux and spectral axis are both specified, check that their lengths
        # match or are off by one (implying the spectral axis stores bin edges)
        if flux is not None and spectral_axis is not None:
            if spectral_axis.shape[0] == flux.shape[-1]:
                if bin_specification == "edges":
                    raise ValueError("A spectral axis input as bin edges"
                        "must have length one greater than the flux axis")
                bin_specification = "centers"
            elif spectral_axis.shape[0] == flux.shape[-1]+1:
                if bin_specification == "centers":
                    raise ValueError("A spectral axis input as bin centers"
                        "must be the same length as the flux axis")
                bin_specification = "edges"
            else:
                raise ValueError(
                    "Spectral axis length ({}) must be the same size or one "
                    "greater (if specifying bin edges) than that of the last "
                    "flux axis ({})".format(spectral_axis.shape[0],
                                            flux.shape[-1]))
    
        # If a WCS is provided, check that the spectral axis is last and reorder
        # the arrays if not
        if wcs is not None and hasattr(wcs, "naxis"):
            if wcs.naxis > 1:
                temp_axes = []
                phys_axes = wcs.world_axis_physical_types
                for i in range(len(phys_axes)):
>                   if phys_axes[i][0:2] == "em":
E                   TypeError: 'NoneType' object is not subscriptable

due to the changes in WCS checks from 9574876; the wcs passed to Spectrum1D by this loader is a minimal instance on a log scale without physical types (wcs.world_axis_physical_types = [None, None]).
Apparently the error has been ignored since, as the remote_data test is allowed to fail.

Very straightforward fix would be to simply not pass wcs (possibly for aspcapStar_loader as well?) on initialising Spectrum1D, as spectral_axis has already been constructed explicitly converting to linear scale, and wcs contains no additional information. But those WCS headers with physical_types None seem to be fairly common.

Update: the failure of manga_rss_loader seems related; here the WCS has
world_axis_physical_types = ['em.wl', None]
so there is a valid wavelength axis, but the comparison attempt raises on the second axis - this should probably be caught as
if phys_axes[i] is not None and phys_axes[i][0:2] == "em"

[dhomeier]

Notably aspcapStar_loader which reads in a virtually identical WCS, but with only a single axis (wcs.world_axis_physical_types = [None]), works since the check is only done for wcs.naxis > 1. So in that case the condition
"Input WCS must have exactly one axis with spectral units"
is not tested at all. Maybe some extra modification in Spectrum1D.__init__ is due then to either

1. accept 0 axes with spectral units for 'LOG-LINEAR' WCS (but only if all axes are dimensionless?) or
2. test always for one axis of type "em" even for wcs.naxis == 1 and never pass log-linear WCS to __init__ (which need an explicitly constructed spectral_axis anyway).

Maybe one option would also be to have Spectrum1D.__init__ itself check for 'LOG' in wcs.ctype[i] and do the conversion internally, but that might fail as those WCS do not usually contain the unit information.

[rosteen]

@dhomeier the broken tests were addressed by #839 , are you satisfied that this can be closed? Or do you want some further consideration of the points you brought up in your second comment?

[dhomeier]

@rosteen thanks for the fix; insofar as the bugs are all fixed (first time all tests are passing in months, I believe!), this can be closed.
Having Spectrum1D.__init__ handle log-linear scales itself would perhaps be a nice feature, but definitely a new feature/enhancement, and thus warrant a different issue (and probably some discussion if it's worth the effort).

[rosteen]

Closing, see #840 for follow-up.