What is the intended and actual result of spectrum.with_spectral_unit?

[dhomeier]

When reviewing the handling of units inside of Spectrum1D for #673 (to ensure things like spectral_axis.to('GHz') work as expected), I noticed the with_spectral_unit method, which also has a test in test_spectral_axis_conversions verifying that it functions, but not really what its function should be.
From the docstring one might expect it to return a Spectrum1D copy with spectral_axis converted to the specified unit, but from the code it becomes quickly clear that to only modifies spectrum.wcs. Or rather tries to:

>>> spectrum = Spectrum1D(flux=[26.0, 52.4] * apu.Jy, spectral_axis=[400, 500] * u.AA)
>>> newspec = spectrum.with_spectral_unit('micron')
>>> newspec.wcs
<WCS(output_frame=SpectralFrame, input_frame=CoordinateFrame, forward_transform=Model: SpectralTabular1D
N_inputs: 1
N_outputs: 1
Parameters: 
  points: (array([0, 1]),)
  lookup_table: [400. 500.] Angstrom
  method: linear
  fill_value: nan
  bounds_error: True)>

for a lookup-table WCS. A "true" (CRVAL/CD1_1-constructed) WCS raises

>>> newspec = wcs1dspec.with_spectral_unit('micron')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/derek/lib/python3.8/site-packages/specutils/spectra/spectrum_mixin.py", line 213, in with_spectral_unit
    new_wcs, new_meta = self._new_spectral_wcs(
  File "/Users/derek/lib/python3.8/site-packages/specutils/spectra/spectrum_mixin.py", line 282, in _new_spectral_wcs
    meta['original_unit'] = self.wcs.unit[0]
AttributeError: 'WCS' object has no attribute 'unit'

But the first use case does not seem to have any effect other than adding a

>>> newspec.meta
OrderedDict([('original_unit', Unit("Angstrom"))])

to the copied spectrum, since _new_spectral_wcs internally is creating the new wcs as
gwcs_from_array(self.spectral_axis)
which discards any modifications to wcs.unit made up to that point.

Any insights on the use case for this?

[nmearl]

Good catch. It seems that this is a victim of the transition to the APE 14 implementation. Before, we actually had wrapper classes around gwcs and fits wcs objects that included extra functionality (like exposed unit attributes), but this is no longer the case. We will have to revisit this implementation to fix it up. Clearly, a failure of the unit test coverage...

The intention is to allow changing the stored spectral axis unit, which is maintained by the wcs (whether gwcs or fits wcs). This was a method to easily alter the wcs information to reflect the desired unit change (i.e. it would also alter the values stored in the wcs to reflect the conversion to the new unit).

[dhomeier]

Good catch. It seems that this is a victim of the transition to the APE 14 implementation. Before, we actually had wrapper classes around gwcs and fits wcs objects that included extra functionality (like exposed unit attributes), but this is no longer the case. We will have to revisit this implementation to fix it up. Clearly, a failure of the unit test coverage...

Looks like changing the call to
gwcs_from_array(self.spectral_axis.to(unit))
might do the job. It does work for the resulting spectrum.spectral_axis as well, which was not clear to me beforehand.

[dhomeier]

Still would need to work out what to do with spectra with WCS without lookup_table and unit, like those loaded with wcs1d-fits. Looks like some more updates to wcs_fits.py are due once #632 is reviewed.

[hamogu]

From the docstring, it clearly seems to indicate that I can change spectral units to e.g. velocity. There is even a rest_value that I can use for that. Yet, it fails miserably:

filename = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'
# The spectrum is in the second HDU of this file.
with fits.open(filename) as f:  
    specdata = f[1].data  
from specutils import Spectrum1D
lamb = 10**specdata['loglam'] * u.AA 
flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') 
spec = Spectrum1D(spectral_axis=lamb, flux=flux) 

specnew = spec.with_spectral_unit(unit=u.km/u.s,
                                  velocity_convention='optical',
                                  rest_value=6365 * u.Angstrom)

Instead of a new spectrum in velocity units relative to the wavelength of Halpha, I get:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [50], line 10
      7 flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1') 
      8 spec = Spectrum1D(spectral_axis=lamb, flux=flux) 
---> 10 specnew = spec.with_spectral_unit(unit=u.km/u.s,
     11                                   velocity_convention='optical',
     12                                   rest_value=6365 * u.Angstrom)

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/specutils/spectra/spectrum_mixin.py:211, in OneDSpectrumMixin.with_spectral_unit(self, unit, velocity_convention, rest_value)
    187 def with_spectral_unit(self, unit, velocity_convention=None,
    188                        rest_value=None):
    189     """
    190     Returns a new spectrum with a different spectral axis unit.
    191 
   (...)
    209 
    210     """
--> 211     new_wcs, new_meta = self._new_spectral_wcs(
    212         unit=unit,
    213         velocity_convention=velocity_convention or self.velocity_convention,
    214         rest_value=rest_value or self.rest_value)
    216     spectrum = self.__class__(flux=self.flux, wcs=new_wcs, meta=new_meta)
    218     return spectrum

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/specutils/spectra/spectrum_mixin.py:274, in OneDSpectrumMixin._new_spectral_wcs(self, unit, velocity_convention, rest_value)
    272 if velocity_convention is not None:
    273     equiv = getattr(u, 'doppler_{0}'.format(velocity_convention))
--> 274     rest_value.to(unit, equivalencies=equiv)
    276 # Store the original unit information for posterity
    277 meta = self._meta.copy()

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/quantity.py:846, in Quantity.to(self, unit, equivalencies, copy)
    842 unit = Unit(unit)
    843 if copy:
    844     # Avoid using to_value to ensure that we make a copy. We also
    845     # don't want to slow down this method (esp. the scalar case).
--> 846     value = self._to_value(unit, equivalencies)
    847 else:
    848     # to_value only copies if necessary
    849     value = self.to_value(unit, equivalencies)

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/quantity.py:800, in Quantity._to_value(self, unit, equivalencies)
    797     equivalencies = self._equivalencies
    798 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    799     # Standard path, let unit to do work.
--> 800     return self.unit.to(unit, self.view(np.ndarray),
    801                         equivalencies=equivalencies)
    803 else:
    804     # The .to() method of a simple unit cannot convert a structured
    805     # dtype, so we work around it, by recursing.
    806     # TODO: deprecate this?
    807     # Convert simple to Structured on initialization?
    808     result = np.empty_like(self.view(np.ndarray))

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/core.py:1130, in UnitBase.to(self, other, value, equivalencies)
   1128     return UNITY
   1129 else:
-> 1130     return self._get_converter(Unit(other),
   1131                                equivalencies=equivalencies)(value)

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/core.py:1047, in UnitBase._get_converter(self, other, equivalencies)
   1044 # if that doesn't work, maybe we can do it with equivalencies?
   1045 try:
   1046     return self._apply_equivalencies(
-> 1047         self, other, self._normalize_equivalencies(equivalencies))
   1048 except UnitsError as exc:
   1049     # Last hope: maybe other knows how to do it?
   1050     # We assume the equivalencies have the unit itself as first item.
   1051     # TODO: maybe better for other to have a `_back_converter` method?
   1052     if hasattr(other, 'equivalencies'):

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/core.py:769, in UnitBase._normalize_equivalencies(equivalencies)
    748 @staticmethod
    749 def _normalize_equivalencies(equivalencies):
    750     """
    751     Normalizes equivalencies, ensuring each is a 4-tuple of the form::
    752 
   (...)
    767     ValueError if an equivalency cannot be interpreted
    768     """
--> 769     normalized = _normalize_equivalencies(equivalencies)
    770     if equivalencies is not None:
    771         normalized += get_current_unit_registry().equivalencies

File ~/mambaforge/envs/kitchensink/lib/python3.10/site-packages/astropy/units/core.py:83, in _normalize_equivalencies(equivalencies)
     79     return []
     81 normalized = []
---> 83 for i, equiv in enumerate(equivalencies):
     84     if len(equiv) == 2:
     85         funit, tunit = equiv

TypeError: 'function' object is not iterable

[rosteen]

Closed by #1119.