From 64393275631b1b4144b2f65d394882ee2be10139 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Mon, 22 Jan 2024 17:18:57 -0500
Subject: [PATCH 01/15] Fix with_spectral_axis when Spectrum1D instantiated
 with a spectral_axis

---
 specutils/spectra/spectrum_mixin.py | 43 ++++++++++++++++++-----------
 specutils/tests/test_spectrum1d.py  |  7 +++--
 specutils/utils/wcs_utils.py        | 22 +++++++++------
 3 files changed, 46 insertions(+), 26 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index e067df957..2d18a6b76 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -5,7 +5,7 @@
 from astropy import units as u
 from astropy.utils.decorators import deprecated
 
-from specutils.utils.wcs_utils import gwcs_from_array
+from specutils.utils.wcs_utils import gwcs_from_array, SpectralGWCS
 
 DOPPLER_CONVENTIONS = {}
 DOPPLER_CONVENTIONS['radio'] = u.doppler_radio
@@ -208,12 +208,31 @@ def with_spectral_unit(self, unit, velocity_convention=None,
                      even if your spectrum has air wavelength units
 
         """
-        new_wcs, new_meta = self._new_spectral_wcs(
-            unit=unit,
-            velocity_convention=velocity_convention or self.velocity_convention,
-            rest_value=rest_value or self.rest_value)
+        velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention # noqa
+        rest_value = rest_value if rest_value is not None else self.rest_value
 
-        spectrum = self.__class__(flux=self.flux, wcs=new_wcs, meta=new_meta)
+        # Store the original unit information for posterity
+        meta = self._meta.copy()
+        orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
+        if 'original_unit' not in self._meta:
+            meta['original_unit'] = orig_unit
+
+        # In this case the original Spectrum1D was instantiated with a spectral_axis and no wcs
+        if isinstance(self.wcs, SpectralGWCS):
+            new_spectral_axis = self.spectral_axis.to(unit,
+                                                      doppler_convention=velocity_convention,
+                                                      doppler_rest=rest_value)
+
+            return self.__class__(flux=self.flux, spectral_axis=new_spectral_axis, meta=meta,
+                                  uncertainty=self.uncertainty, mask=self.mask)
+
+        # Otherwise we handle the original wcs properly
+        new_wcs = self._new_spectral_wcs(unit=unit, orig_unit=orig_unit,
+                                         velocity_convention=velocity_convention,
+                                         rest_value=rest_value)
+
+        spectrum = self.__class__(flux=self.flux, wcs=new_wcs, meta=meta,
+                                  uncertainty=self.uncertainty, mask=self.mask)
 
         return spectrum
 
@@ -242,7 +261,7 @@ def _new_wcs_argument_validation(self, unit, velocity_convention,
 
         return unit
 
-    def _new_spectral_wcs(self, unit, velocity_convention=None,
+    def _new_spectral_wcs(self, unit, orig_unit, velocity_convention=None,
                           rest_value=None):
         """
         Returns a new WCS with a different Spectral Axis unit.
@@ -273,18 +292,10 @@ def _new_spectral_wcs(self, unit, velocity_convention=None,
             equiv = getattr(u, 'doppler_{0}'.format(velocity_convention))
             rest_value.to(unit, equivalencies=equiv)
 
-        # Store the original unit information for posterity
-        meta = self._meta.copy()
-
-        orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
-
-        if 'original_unit' not in self._meta:
-            meta['original_unit'] = orig_unit
-
         # Create the new wcs object
         if isinstance(unit, u.UnitBase) and unit.is_equivalent(
                 orig_unit, equivalencies=u.spectral()):
-            return gwcs_from_array(self.spectral_axis), meta
+            return gwcs_from_array(self.spectral_axis)
 
         raise u.UnitConversionError(f"WCS units incompatible: {unit} and {orig_unit}.")
 
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 07efa7fbb..81727eb86 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -150,10 +150,13 @@ def test_spectral_axis_conversions():
     with pytest.raises(ValueError):
         spec.velocity
 
-    spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm,
+    spec = Spectrum1D(spectral_axis=np.arange(100, 150) * u.nm,
                       flux=np.random.randn(49) * u.Jy)
 
-    new_spec = spec.with_spectral_unit(u.GHz)  # noqa
+    new_spec = spec.with_spectral_unit(u.km/u.s, rest_value=125*u.um,
+                                       velocity_convention="relativistic")
+
+    assert new_spec.spectral_axis.unit == u.km/u.s
 
 
 def test_spectral_slice():
diff --git a/specutils/utils/wcs_utils.py b/specutils/utils/wcs_utils.py
index 24b315758..dce020781 100644
--- a/specutils/utils/wcs_utils.py
+++ b/specutils/utils/wcs_utils.py
@@ -8,6 +8,18 @@
 from gwcs import coordinate_frames as cf
 
 
+class SpectralGWCS(GWCS):
+    def __init__(self, *args, **kwargs):
+        self.original_unit = kwargs.pop("original_unit")
+        super().__init__(*args, **kwargs)
+
+    def pixel_to_world(self, *args, **kwargs):
+        if self.original_unit == '':
+            return u.Quantity(super().pixel_to_world_values(*args, **kwargs))
+        return super().pixel_to_world(*args, **kwargs).to(
+            self.original_unit, equivalencies=u.spectral())
+
+
 def refraction_index(wavelength, method='Griesen2006', co2=None):
     """
     Calculates the index of refraction of dry air at standard temperature
@@ -210,14 +222,8 @@ def gwcs_from_array(array):
         forward_transform.inverse = SpectralTabular1D(
                 array[::-1], lookup_table=np.arange(len(array))[::-1])
 
-    class SpectralGWCS(GWCS):
-        def pixel_to_world(self, *args, **kwargs):
-            if orig_array.unit == '':
-                return u.Quantity(super().pixel_to_world_values(*args, **kwargs))
-            return super().pixel_to_world(*args, **kwargs).to(
-                orig_array.unit, equivalencies=u.spectral())
-
-    tabular_gwcs = SpectralGWCS(forward_transform=forward_transform,
+    tabular_gwcs = SpectralGWCS(original_unit = orig_array.unit,
+                                forward_transform=forward_transform,
                                 input_frame=coord_frame,
                                 output_frame=spec_frame)
 

From b8f439c78994118fd1fd49980ec3813e92ae7116 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Mon, 22 Jan 2024 18:42:44 -0500
Subject: [PATCH 02/15] Deprecate old method name, new one is more specific

---
 specutils/spectra/spectrum_mixin.py | 6 ++++++
 specutils/tests/test_spectrum1d.py  | 2 +-
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 2d18a6b76..8f6ea461d 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -184,8 +184,14 @@ def velocity(self):
 
         return new_data
 
+    @deprecated('v1.13', alternative="with_spectral_axis_unit")
     def with_spectral_unit(self, unit, velocity_convention=None,
                            rest_value=None):
+        self.with_spectral_axis_unit(unit, velocity_convention=velocity_convention,
+                                     rest_value=rest_value)
+
+    def with_spectral_axis_unit(self, unit, velocity_convention=None,
+                           rest_value=None):
         """
         Returns a new spectrum with a different spectral axis unit.
 
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index 81727eb86..e3c35abd8 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -153,7 +153,7 @@ def test_spectral_axis_conversions():
     spec = Spectrum1D(spectral_axis=np.arange(100, 150) * u.nm,
                       flux=np.random.randn(49) * u.Jy)
 
-    new_spec = spec.with_spectral_unit(u.km/u.s, rest_value=125*u.um,
+    new_spec = spec.with_spectral_axis_unit(u.km/u.s, rest_value=125*u.um,
                                        velocity_convention="relativistic")
 
     assert new_spec.spectral_axis.unit == u.km/u.s

From f68ab6e8ef701922d3728459a94a08db3dcd841c Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Mon, 22 Jan 2024 19:07:02 -0500
Subject: [PATCH 03/15] Works for WCS-initialized spectra as well

---
 specutils/spectra/spectrum_mixin.py | 73 +++++++----------------------
 1 file changed, 16 insertions(+), 57 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 8f6ea461d..032d4c709 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -193,7 +193,10 @@ def with_spectral_unit(self, unit, velocity_convention=None,
     def with_spectral_axis_unit(self, unit, velocity_convention=None,
                            rest_value=None):
         """
-        Returns a new spectrum with a different spectral axis unit.
+        Returns a new spectrum with a different spectral axis unit. Note that this creates a new
+        object using the converted spectral axis and thus drops the original WCS, if it existed.
+        The original WCS will be stored in the ``original_wcs`` entry of the new object's ``meta``
+        dictionary.
 
         Parameters
         ----------
@@ -214,33 +217,27 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
                      even if your spectrum has air wavelength units
 
         """
+
         velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention # noqa
         rest_value = rest_value if rest_value is not None else self.rest_value
+        unit = self._new_wcs_argument_validation(unit, velocity_convention,
+                                                 rest_value)
 
-        # Store the original unit information for posterity
+        # Store the original unit information and WCS for posterity
         meta = self._meta.copy()
-        orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
+
         if 'original_unit' not in self._meta:
+            orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
             meta['original_unit'] = orig_unit
 
-        # In this case the original Spectrum1D was instantiated with a spectral_axis and no wcs
-        if isinstance(self.wcs, SpectralGWCS):
-            new_spectral_axis = self.spectral_axis.to(unit,
-                                                      doppler_convention=velocity_convention,
-                                                      doppler_rest=rest_value)
+        if 'original_wcs' not in self.meta:
+            meta['original_wcs'] = self.wcs.copy()
 
-            return self.__class__(flux=self.flux, spectral_axis=new_spectral_axis, meta=meta,
-                                  uncertainty=self.uncertainty, mask=self.mask)
+        new_spectral_axis = self.spectral_axis.to(unit, doppler_convention=velocity_convention,
+                                                  doppler_rest=rest_value)
 
-        # Otherwise we handle the original wcs properly
-        new_wcs = self._new_spectral_wcs(unit=unit, orig_unit=orig_unit,
-                                         velocity_convention=velocity_convention,
-                                         rest_value=rest_value)
-
-        spectrum = self.__class__(flux=self.flux, wcs=new_wcs, meta=meta,
-                                  uncertainty=self.uncertainty, mask=self.mask)
-
-        return spectrum
+        return self.__class__(flux=self.flux, spectral_axis=new_spectral_axis, meta=meta,
+                              uncertainty=self.uncertainty, mask=self.mask)
 
     def _new_wcs_argument_validation(self, unit, velocity_convention,
                                      rest_value):
@@ -267,44 +264,6 @@ def _new_wcs_argument_validation(self, unit, velocity_convention,
 
         return unit
 
-    def _new_spectral_wcs(self, unit, orig_unit, velocity_convention=None,
-                          rest_value=None):
-        """
-        Returns a new WCS with a different Spectral Axis unit.
-
-        Parameters
-        ----------
-        unit : :class:`~astropy.units.Unit`
-            Any valid spectral unit: velocity, (wave)length, or frequency.
-            Only vacuum units are supported.
-        velocity_convention : 'relativistic', 'radio', or 'optical'
-            The velocity convention to use for the output velocity axis.
-            Required if the output type is velocity. This can be either one
-            of the above strings, or an `astropy.units` equivalency.
-        rest_value : :class:`~astropy.units.Quantity`
-            A rest wavelength or frequency with appropriate units.  Required if
-            output type is velocity.  The cube's WCS should include this
-            already if the *input* type is velocity, but the WCS's rest
-            wavelength/frequency can be overridden with this parameter.
-
-            .. note: This must be the rest frequency/wavelength *in vacuum*,
-                     even if your cube has air wavelength units
-        """
-
-        unit = self._new_wcs_argument_validation(unit, velocity_convention,
-                                                 rest_value)
-
-        if velocity_convention is not None:
-            equiv = getattr(u, 'doppler_{0}'.format(velocity_convention))
-            rest_value.to(unit, equivalencies=equiv)
-
-        # Create the new wcs object
-        if isinstance(unit, u.UnitBase) and unit.is_equivalent(
-                orig_unit, equivalencies=u.spectral()):
-            return gwcs_from_array(self.spectral_axis)
-
-        raise u.UnitConversionError(f"WCS units incompatible: {unit} and {orig_unit}.")
-
     def _check_strictly_increasing_decreasing(self):
         """
         Check that the self._spectral_axis is strictly increasing or decreasing

From 86f474458e464f23829ae0c1db1c792fbf5583ad Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Mon, 22 Jan 2024 19:09:15 -0500
Subject: [PATCH 04/15] Remove unneeded imports

---
 specutils/spectra/spectrum_mixin.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 032d4c709..be2b740fd 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -5,8 +5,6 @@
 from astropy import units as u
 from astropy.utils.decorators import deprecated
 
-from specutils.utils.wcs_utils import gwcs_from_array, SpectralGWCS
-
 DOPPLER_CONVENTIONS = {}
 DOPPLER_CONVENTIONS['radio'] = u.doppler_radio
 DOPPLER_CONVENTIONS['optical'] = u.doppler_optical

From 88b2abd9409083a2128bb8204c8d6866ee396176 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Mon, 22 Jan 2024 19:12:02 -0500
Subject: [PATCH 05/15] Add kwarg default

---
 specutils/utils/wcs_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/utils/wcs_utils.py b/specutils/utils/wcs_utils.py
index dce020781..6a0ec8e96 100644
--- a/specutils/utils/wcs_utils.py
+++ b/specutils/utils/wcs_utils.py
@@ -10,7 +10,7 @@
 
 class SpectralGWCS(GWCS):
     def __init__(self, *args, **kwargs):
-        self.original_unit = kwargs.pop("original_unit")
+        self.original_unit = kwargs.pop("original_unit", "")
         super().__init__(*args, **kwargs)
 
     def pixel_to_world(self, *args, **kwargs):

From 0ff7bac09dd5ad7efbe17816c48e825cb7549b4d Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Tue, 23 Jan 2024 09:46:03 -0500
Subject: [PATCH 06/15] Add copy and deepcopy methods

Codestyle
---
 specutils/utils/wcs_utils.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/specutils/utils/wcs_utils.py b/specutils/utils/wcs_utils.py
index 6a0ec8e96..cce7fb609 100644
--- a/specutils/utils/wcs_utils.py
+++ b/specutils/utils/wcs_utils.py
@@ -13,6 +13,28 @@ def __init__(self, *args, **kwargs):
         self.original_unit = kwargs.pop("original_unit", "")
         super().__init__(*args, **kwargs)
 
+    def copy(self):
+        """
+        Return a shallow copy of the object.
+
+        Convenience method so user doesn't have to import the
+        :mod:`copy` stdlib module.
+
+        .. warning::
+            Use `deepcopy` instead of `copy` unless you know why you need a
+            shallow copy.
+        """
+        return copy.copy(self)
+
+    def deepcopy(self):
+        """
+        Return a deep copy of the object.
+
+        Convenience method so user doesn't have to import the
+        :mod:`copy` stdlib module.
+        """
+        return copy.deepcopy(self)
+
     def pixel_to_world(self, *args, **kwargs):
         if self.original_unit == '':
             return u.Quantity(super().pixel_to_world_values(*args, **kwargs))

From 5e9cc7f356ec94a8c2168931bad06a1081dca55f Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>
Date: Tue, 23 Jan 2024 11:22:56 -0500
Subject: [PATCH 07/15] Update specutils/spectra/spectrum_mixin.py

Co-authored-by: Derek Homeier <709020+dhomeier@users.noreply.github.com>
---
 specutils/spectra/spectrum_mixin.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index be2b740fd..02d2c3fe7 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -229,7 +229,7 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
             meta['original_unit'] = orig_unit
 
         if 'original_wcs' not in self.meta:
-            meta['original_wcs'] = self.wcs.copy()
+            meta['original_wcs'] = self.wcs.deepcopy()
 
         new_spectral_axis = self.spectral_axis.to(unit, doppler_convention=velocity_convention,
                                                   doppler_rest=rest_value)

From 5221ce4c60657a976b7623cd36055fe49b5be902 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Tue, 23 Jan 2024 11:25:16 -0500
Subject: [PATCH 08/15] Add changelog entry

---
 CHANGES.rst | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 76a1098af..15ba630c4 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -10,6 +10,9 @@ Bug Fixes
 - ``template_correlate`` no longer errors when used on a ``Spectrum1D`` that lacks an
   ``uncertainty`` array. [#1118]
 
+- ``with_spectral_unit`` has been changed to ``with_spectral_axis_unit`` and actually works
+  now. [#1119]
+
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -64,7 +67,7 @@ New Features
 
 - ``wcs1d-fits`` loader now reads and writes celestial components of
   of multi-dimensional WCS, and handles ``mask`` and ``uncertainty``
-  attributes. [#1009] 
+  attributes. [#1009]
 
 - Added support for reading from files with flux in counts. [#1018]
 
@@ -183,7 +186,7 @@ Bug Fixes
   ``Spectrum1D`` without WCS nor spectral axis and the spatial-spatial dimension
   is smaller than spectral dimension. [#926]
 
-- Fixed WCS not accurately reflecting the updated spectral axis after slicing a 
+- Fixed WCS not accurately reflecting the updated spectral axis after slicing a
   ``Spectrum1D``. [#918]
 
 Other Changes and Additions
@@ -212,8 +215,8 @@ New Features
 - Convolution-based smoothing will now apply a 1D kernel to multi-dimensional fluxes
   by convolving along the spectral axis only, rather than raising an error. [#885]
 
-- ``template_comparison`` now handles ``astropy.nddata.Variance`` and 
-  ``astropy.nddata.InverseVariance`` uncertainties instead of assuming 
+- ``template_comparison`` now handles ``astropy.nddata.Variance`` and
+  ``astropy.nddata.InverseVariance`` uncertainties instead of assuming
   the uncertainty is standard deviation. [#899]
 
 Bug Fixes
@@ -253,7 +256,7 @@ New Features
 
 - Allow overriding existing astropy registry elements. [#861]
 
-- ``Spectrum1D`` will now swap the spectral axis with the last axis on initialization 
+- ``Spectrum1D`` will now swap the spectral axis with the last axis on initialization
   if it can be identified from the WCS and is not last, rather than erroring. [#654, #822]
 
 Bug Fixes

From c54a30ee08b826ed528acf6b1c4d3735dea67e65 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Tue, 23 Jan 2024 11:30:03 -0500
Subject: [PATCH 09/15] Test to make sure we kept the old WCS in meta

---
 specutils/tests/test_spectrum1d.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index e3c35abd8..f9a8db12c 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -157,6 +157,9 @@ def test_spectral_axis_conversions():
                                        velocity_convention="relativistic")
 
     assert new_spec.spectral_axis.unit == u.km/u.s
+    assert new_spec.wcs.world_axis_units[0] == "km.s**-1"
+    # Make sure meta stored the old WCS correctly
+    assert new_spec.meta["original_wcs"].world_axis_units[0] == "nm"
 
 
 def test_spectral_slice():

From c8159d25253ba2757f8011d60151bd01200a88c8 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Tue, 23 Jan 2024 11:31:51 -0500
Subject: [PATCH 10/15] Better name for meta attribute, add test for it

---
 specutils/spectra/spectrum_mixin.py | 4 ++--
 specutils/tests/test_spectrum1d.py  | 1 +
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 02d2c3fe7..85333f923 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -224,9 +224,9 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
         # Store the original unit information and WCS for posterity
         meta = self._meta.copy()
 
-        if 'original_unit' not in self._meta:
+        if 'original_spectral_axis_unit' not in self._meta:
             orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
-            meta['original_unit'] = orig_unit
+            meta['original_spectral_axis_unit'] = orig_unit
 
         if 'original_wcs' not in self.meta:
             meta['original_wcs'] = self.wcs.deepcopy()
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index f9a8db12c..ba9060ac1 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -160,6 +160,7 @@ def test_spectral_axis_conversions():
     assert new_spec.wcs.world_axis_units[0] == "km.s**-1"
     # Make sure meta stored the old WCS correctly
     assert new_spec.meta["original_wcs"].world_axis_units[0] == "nm"
+    assert new_spec.meta["original_spectral_axis_unit"] == "nm"
 
 
 def test_spectral_slice():

From 4949bea8660fbdb52deba7e3b3e135ecbc30186e Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>
Date: Tue, 30 Jan 2024 12:59:29 -0500
Subject: [PATCH 11/15] Apply suggestions from code review

Co-authored-by: Derek Homeier <709020+dhomeier@users.noreply.github.com>
---
 specutils/spectra/spectrum_mixin.py |  2 +-
 specutils/tests/test_spectrum1d.py  | 13 +++++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 85333f923..5a358bb26 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -222,7 +222,7 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
                                                  rest_value)
 
         # Store the original unit information and WCS for posterity
-        meta = self._meta.copy()
+        meta = deepcopy(self._meta)
 
         if 'original_spectral_axis_unit' not in self._meta:
             orig_unit = self.wcs.unit[0] if hasattr(self.wcs, 'unit') else self.spectral_axis.unit
diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index ba9060ac1..cdf342ae5 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -162,6 +162,19 @@ def test_spectral_axis_conversions():
     assert new_spec.meta["original_wcs"].world_axis_units[0] == "nm"
     assert new_spec.meta["original_spectral_axis_unit"] == "nm"
 
+    wcs_dict = {"CTYPE1": "WAVE", "CRVAL1": 3.622e3, "CDELT1": 8e-2,
+                "CRPIX1": 0, "CUNIT1": "Angstrom"}
+    wcs_spec = Spectrum1D(flux=np.random.randn(49) * u.Jy, wcs=WCS(wcs_dict),
+                          meta={'header': wcs_dict.copy()})
+    new_spec = wcs_spec.with_spectral_axis_unit(u.km/u.s, rest_value=125*u.um,
+                                       velocity_convention="relativistic")
+    new_spec.meta['original_wcs'].wcs.crval = [3.777e-7]
+    new_spec.meta['header']['CRVAL1'] = 3777.0
+
+    assert wcs_spec.wcs.wcs.crval[0] == 3.622e-7
+    assert wcs_spec.meta['header']['CRVAL1'] == 3622.
+
+
 
 def test_spectral_slice():
     spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm,

From 2a4f39aaaffe9993affe4a1131ed0c9021f2453c Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Tue, 30 Jan 2024 13:03:26 -0500
Subject: [PATCH 12/15] Fix codestyle

---
 specutils/tests/test_spectrum1d.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py
index cdf342ae5..eea2a43ad 100644
--- a/specutils/tests/test_spectrum1d.py
+++ b/specutils/tests/test_spectrum1d.py
@@ -175,7 +175,6 @@ def test_spectral_axis_conversions():
     assert wcs_spec.meta['header']['CRVAL1'] == 3622.
 
 
-
 def test_spectral_slice():
     spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm,
                       flux=np.random.random(10) * u.Jy)

From eea7796137af22b223e0ddcab9dde426bd5c268c Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>
Date: Thu, 1 Feb 2024 13:04:22 -0500
Subject: [PATCH 13/15] Update docstring

Co-authored-by: Derek Homeier <709020+dhomeier@users.noreply.github.com>
---
 specutils/spectra/spectrum_mixin.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 5a358bb26..717cae650 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -192,9 +192,9 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
                            rest_value=None):
         """
         Returns a new spectrum with a different spectral axis unit. Note that this creates a new
-        object using the converted spectral axis and thus drops the original WCS, if it existed.
-        The original WCS will be stored in the ``original_wcs`` entry of the new object's ``meta``
-        dictionary.
+        object using the converted spectral axis and thus drops the original WCS, if it existed,
+        replacing it with a :class:`~specutils.utils.wcs_utils.SpectralGWCS`. The original WCS
+        will be stored in the ``original_wcs`` entry of the new object's ``meta`` dictionary.
 
         Parameters
         ----------

From 217dbf0bf765e571cc38fc6785eb799e7315f1ce Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Thu, 1 Feb 2024 15:16:22 -0500
Subject: [PATCH 14/15] Change docstring to fix docs build failure. Also fixed
 some incorrect references to SpectralRegion

---
 specutils/analysis/uncertainty.py   | 8 ++++----
 specutils/spectra/spectrum_mixin.py | 5 +++--
 specutils/utils/wcs_utils.py        | 4 ++++
 3 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/specutils/analysis/uncertainty.py b/specutils/analysis/uncertainty.py
index edd6a99a3..2ee6b110c 100644
--- a/specutils/analysis/uncertainty.py
+++ b/specutils/analysis/uncertainty.py
@@ -22,7 +22,7 @@ def snr(spectrum, region=None):
     spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
         The spectrum object overwhich the equivalent width will be calculated.
 
-    region: `~specutils.utils.SpectralRegion` or list of `~specutils.utils.SpectralRegion`
+    region: `~specutils.SpectralRegion` or list of `~specutils.SpectralRegion`
         Region within the spectrum to calculate the SNR.
 
     Returns
@@ -67,7 +67,7 @@ def _snr_single_region(spectrum, region=None):
     spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
         The spectrum object overwhich the equivalent width will be calculated.
 
-    region: `~specutils.utils.SpectralRegion`
+    region: `~specutils.SpectralRegion`
         Region within the spectrum to calculate the SNR.
 
     Returns
@@ -109,7 +109,7 @@ def snr_derived(spectrum, region=None):
     spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
         The spectrum object overwhich the equivalent width will be calculated.
 
-    region: `~specutils.utils.SpectralRegion`
+    region: `~specutils.SpectralRegion`
         Region within the spectrum to calculate the SNR.
 
     Returns
@@ -155,7 +155,7 @@ def _snr_derived(spectrum, region=None):
     spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
         The spectrum object overwhich the equivalent width will be calculated.
 
-    region: `~specutils.utils.SpectralRegion`
+    region: `~specutils.SpectralRegion`
         Region within the spectrum to calculate the SNR.
 
     Returns
diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py
index 717cae650..f57a5c595 100644
--- a/specutils/spectra/spectrum_mixin.py
+++ b/specutils/spectra/spectrum_mixin.py
@@ -193,8 +193,9 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None,
         """
         Returns a new spectrum with a different spectral axis unit. Note that this creates a new
         object using the converted spectral axis and thus drops the original WCS, if it existed,
-        replacing it with a :class:`~specutils.utils.wcs_utils.SpectralGWCS`. The original WCS
-        will be stored in the ``original_wcs`` entry of the new object's ``meta`` dictionary.
+        replacing it with a lookup-table :class:`~gwcs.wcs.WCS` based on the new spectral axis. The
+        original WCS will be stored in the ``original_wcs`` entry of the new object's ``meta``
+        dictionary.
 
         Parameters
         ----------
diff --git a/specutils/utils/wcs_utils.py b/specutils/utils/wcs_utils.py
index cce7fb609..bb4a5f38e 100644
--- a/specutils/utils/wcs_utils.py
+++ b/specutils/utils/wcs_utils.py
@@ -9,6 +9,10 @@
 
 
 class SpectralGWCS(GWCS):
+    """
+    This is a placeholder lookup-table GWCS created when a :class:`~specutils.Spectrum1D` is
+    instantiated with a ``spectral_axis`` and no WCS.
+    """
     def __init__(self, *args, **kwargs):
         self.original_unit = kwargs.pop("original_unit", "")
         super().__init__(*args, **kwargs)

From 2e357c8649ca5a4d0d43dbb10b1853dbc08b05e7 Mon Sep 17 00:00:00 2001
From: Ricky O'Steen <rosteen@stsci.edu>
Date: Thu, 1 Feb 2024 15:23:35 -0500
Subject: [PATCH 15/15] Fix one more SpectralRegion reference

---
 specutils/analysis/moment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specutils/analysis/moment.py b/specutils/analysis/moment.py
index 224ae4b31..4b50cbb5d 100644
--- a/specutils/analysis/moment.py
+++ b/specutils/analysis/moment.py
@@ -21,7 +21,7 @@ def moment(spectrum, regions=None, order=0, axis=-1):
     spectrum : `~specutils.spectra.spectrum1d.Spectrum1D`
         The spectrum object over which the width will be calculated.
 
-    regions: `~specutils.utils.SpectralRegion` or list of `~specutils.utils.SpectralRegion`
+    regions: `~specutils.SpectralRegion` or list of `~specutils.SpectralRegion`
         Region within the spectrum to calculate the gaussian sigma width. If
         regions is `None`, computation is performed over entire spectrum.