Unable to write Spectrum1D as wcs1d-fits after using shift_spectrum_to()

[bhgr216]

After shifting a spectrum, attempting to write a Spectrum1D to a fits file with the format wcs1d-fits using Spectrum1D.write(format="wcs1d-fits") raises the following error: IndexError: index 6 is out of bounds for axis 0 with size 1. Using the same function with format="tabular-fits" works as expected.

I suspect this is because the wcs is not updated after the shift:

print(spectrum.wcs)
spectrum.shift_spectrum_to(redshift=0.5)
print(spectrum.wcs)

Output:

WCS Keywords

Number of WCS axes: 1
CTYPE : 'LINEAR' 
CRVAL : 5621.3397216801 
CRPIX : 1.0 
PC1_1  : 1.0 
CDELT : 1.577376238474 
NAXIS : 0  0
WCS Keywords

Number of WCS axes: 1
CTYPE : 'LINEAR' 
CRVAL : 5621.3397216801 
CRPIX : 1.0 
PC1_1  : 1.0 
CDELT : 1.577376238474 
NAXIS : 0  0

[bhgr216]

I now realize this is because after correcting for redshift/radial velocity the spectral axis will not have a constant delta, which doesn't work with a fits header. A solution would be to allow wcs to recognize radial velocity as a header keyword, which would allow us to at least store the value in the fits header while keeping the original spectral axis, but that's an issue for AstroPy rather than specutils.

Alternatively, Spectrum1D.write() could add a RADVEL keyword to the header manually after using AstroPy to create the fits file (which I assume is how it works).