redshift applied in the wrong direction in Spectrum1D object

[stephjuneau]

Hi, I just tried filling in the redshift value for a Spectrum1D object after the redshift information became available, and the spectral axis was immediately multiplied by (1+z). I would normally expect that the wavelength will be reported in observed frame by default and not changed when the user adds the measured redshift to the object but that the redshift allows the user to convert from observed frame to rest-frame. As a user, I would also expect that if there is a "to_rest" operation, it will apply to the spectral axis and to the flux and uncertainty as well, taking account the units to treat the f_lambda and f_nu units properly (the latter is already in issue #741 #741). I would also hope to be able to know whether the wavelength/fluxes are in rest or observed frame (maybe a keyword with rest=True or False ?).

Here is a simple code that I used in a Notebook based on a demo ran by Erik Tollerud:

from astropy.utils.data import download_file
from specutils import Spectrum1D
%matplotlib inline
import matplotlib.pyplot as plt

spec_fn = download_file('https://stsci.box.com/shared/static/b22b1fzhimtdqfp8597m4bg67kovvauu.fits', cache=True)

jwst_spec = Spectrum1D.read(spec_fn)
plt.plot(jwst_spec.spectral_axis, jwst_spec.flux)

jwst_spec.redshift = 1.5
plt.plot(jwst_spec.spectral_axis, jwst_spec.flux, color='r')

[eteq]

Thanks @stephjuneau ! I think I worked out a more minimal example that demonstrates the issue :

>>> wl = np.linspace(5000, 9000, 100)*u.angstrom
>>> s = specutils.Spectrum1D(spectral_axis=wl, flux=np.ones(100)*u.Jy)
>>> s.spectral_axis
[5000, ..., 9000]A˚
>>> s.redshift = 0.5
>>> s.spectral_axis.to_rest()
[5000, ..., 9000]A˚

Whereas @stephjuneau was expecting this to yield a rest wavelength from 3333 to 6000.

however, this is also true:

>>> s = specutils.Spectrum1D(spectral_axis=wl, flux=np.ones(100)*u.Jy)
>>> s.redshift = 0.5
>>> s.spectral_axis
[7500,  ..., 13500]A˚

Which is the behavior @stephjuneau says above is happening.

But to get the behavior I think @stephjuneau wants the following works:

>>> s = specutils.Spectrum1D(spectral_axis=x,flux=np.ones(100)*u.Jy, redshift=.5)
>>> s.spectral_axis.to_rest()
[3333.333,  ..., 6000]A˚

that is, if you give the redshift at spectrum creation time it does the right thing. But @stephjuneau is talking about a workflow where you load the spectrum, look at it, and only after that want to assign a redshift.

I think there's three possible solutions to this:

1. Add an additional method set_redshift_to or something like that, which sets the redshift without actually shifting the spectral_axis
2. Change Spectrum1D.redshift to do the above, and add a shift_spectrum_to method that does the current behavior
3. Have an option of Spectrum1D.redshift = None which means "don't know the redshift. Only assign it when I set this the first time. From then on do the shifting of the wavelength axis"

I slightly favor 3 but can be convinced

cc @astrofrog and @nmearl . I feel like we've talked about this before, but I couldn't find an issue about it so...

[eteq]

(Oh, and I know SpectralCoord is implemented in astropy, but I think the behavior I'm describing above makes more sense here in specutils since Spectrum1D actually is a spectrum that users actively want to fiddle with, whereas SpectralCoord is more like an immutable coordinate where once you've set the redshift it's set)