Merge pull request #437 from GreenBankObservatory/401-changing-the-fr…

…ame-of-a-spectrum-to-itself-changes-its-spectral-axis

401 changing the frame of a spectrum to itself changes its spectral axis

notebooks/examples/flagging.ipynb

@@ -514,7 +514,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.0"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,

src/dysh/coordinates/core.py

@@ -54,8 +54,22 @@
     "galactocentric": coord.Galactocentric,
     "topocentric": coord.ITRS,  # but need to add observatory position
     "topo": coord.ITRS,  # but need to add observatory position
+    "itrs": coord.ITRS,  # but need to add observatory positionsc.sp
 }
 
+astropy_convenience_frame_names = {
+    "bary": "icrs",
+    "barycentric": "icrs",
+    "heliocentric": "icrs",
+    "helio": "icrs",
+    "geo": "icrs",
+    "geocentric": "icrs",
+    "topocentric": "itrs",
+    "topo": "itrs",
+    "vlsr": "lsrk",
+}
+
+
 # astropy-sanctioned coordinate frame string to label
 frame_to_label = {
     "itrs": "Topocentric",
@@ -450,7 +464,6 @@ def veltofreq(velocity, restfreq, veldef):
     vdef = veldef_to_convention(veldef)
     if vdef is None:
         raise ValueError(f"Unrecognized VELDEF: {veldef}")
-    vc = velocity / const.c
     if vdef == "radio":
         radio = u.doppler_radio(restfreq)
         frequency = velocity.to(u.Hz, equivalencies=radio)
@@ -469,7 +482,6 @@ def veltofreq(velocity, restfreq, veldef):
 def change_ctype(ctype, toframe):
     # when changing frame, we should also change CTYPE1. Pretty sure GBTIDL does not do this
     prefix = ctype[0:4]
-    postfix = ctype[4:]
     newpostfix = reverse_frame_dict[toframe]
     newctype = prefix + newpostfix
     # print(f"changing {ctype} to {newctype}")

src/dysh/spectra/spectrum.py

@@ -8,7 +8,7 @@
 import astropy.units as u
 import numpy as np
 import pandas as pd
-from astropy.coordinates import SkyCoord, SpectralCoord, StokesCoord
+from astropy.coordinates import ITRS, SkyCoord, SpectralCoord, StokesCoord
 from astropy.coordinates.spectral_coordinate import attach_zero_velocities
 from astropy.io import registry
 from astropy.io.fits.verify import VerifyWarning
@@ -27,6 +27,7 @@
 from ..coordinates import (  # is_topocentric,; topocentric_velocity_to_frame,
     KMS,
     Observatory,
+    astropy_convenience_frame_names,
     astropy_frame_dict,
     change_ctype,
     get_velocity_in_frame,
@@ -35,7 +36,7 @@
     sanitize_skycoord,
     veldef_to_convention,
 )
-from ..log import HistoricalBase, log_call_to_history, logger
+from ..log import HistoricalBase, log_call_to_history  # , logger
 from ..plot import specplot as sp
 from ..util import minimum_string_match
 from . import baseline, get_spectral_equivalency
@@ -722,11 +723,13 @@ def velocity_axis_to(self, unit=KMS, toframe=None, doppler_convention=None):
             The converted spectral axis velocity
         """
         if toframe is not None and toframe != self.velocity_frame:
-            self.set_frame(toframe)
+            s = self.with_frame(toframe)
+        else:
+            s = self
         if doppler_convention is not None:
-            return self._spectral_axis.to(unit=unit, doppler_convention=doppler_convention).to(unit)
+            return s._spectral_axis.to(unit=unit, doppler_convention=doppler_convention).to(unit)
         else:
-            return self.axis_velocity(unit)
+            return s.axis_velocity(unit)
 
     def get_velocity_shift_to(self, toframe):
         if self._target is None:
@@ -735,37 +738,51 @@ def get_velocity_shift_to(self, toframe):
 
     @log_call_to_history
     def set_frame(self, toframe):
-        # @todo VELDEF should be changed as well?
         """Set the sky coordinate and doppler tracking reference frame of this Spectrum. The header 'CTYPE1' will be changed accordingly.
 
         To make a copy of this Spectrum with new coordinate referece frmae instead, use `with_frame`.
 
         Parameters
         ----------
-        toframe - str
-            The coordinate reference frame identifying string, as used by astropy, e.g. 'hcrs', 'icrs', etc.
+        toframe - str, or ~astropy.coordinates.BaseCoordinateFrame, or ~astropy.coordinates.SkyCoord
+            The coordinate reference frame identifying string, as used by astropy, e.g. 'hcrs', 'icrs', etc.,
+            or an actual coordinate system instance
         """
-        if "topo" in toframe:
-            actualframe = self.observer
-        else:
-            actualframe = astropy_frame_dict.get(toframe, toframe)
-        self._spectral_axis = self._spectral_axis.with_observer_stationary_relative_to(actualframe)
+
+        tfl = toframe
+        if isinstance(toframe, str):
+            tfl = toframe.lower()
+            tfl = astropy_convenience_frame_names.get(tfl, tfl)
+            if "itrs" in tfl:
+                if isinstance(self._observer, ITRS):
+                    return  # nothing to be done, we already have the correct axis
+                raise ValueError(
+                    "For topographic or ITRS coordaintes, you must supply a full astropy Coordinate instance."
+                )
+            elif self._velocity_frame == tfl:
+                return  # the frame is already the requested frame
+
+        self._spectral_axis = self._spectral_axis.with_observer_stationary_relative_to(tfl)
+        self._observer = self._spectral_axis.observer
         # This line is commented because:
         # SDFITS defines CTYPE1 as always being the TOPO frequency.
         # See Issue #373 on GitHub.
         # self._meta["CTYPE1"] = change_ctype(self._meta["CTYPE1"], toframe)
-        if isinstance(actualframe, str):
-            self._velocity_frame = actualframe
+        if isinstance(tfl, str):
+            self._velocity_frame = tfl
         else:
-            self._velocity_frame = actualframe.name
+            self._velocity_frame = tfl.name
+        # While it is incorrect to change CTYPE1, it is reasonable to change VELDEF
+        self.meta["VELDEF"] = change_ctype(self.meta["VELDEF"], self._velocity_frame)
 
     def with_frame(self, toframe):
         """Return a copy of this Spectrum with a new coordinate reference frame.
 
         Parameters
         ----------
-        toframe - str
-            The coordinate reference frame identifying string, as used by astropy, e.g. 'hcrs', 'icrs', etc.
+        toframe - str, ~astropy.coordinates.BaseCoordinateFrame, or ~astropy.coordinates.SkyCoord
+            The coordinate reference frame identifying string, as used by astropy, e.g. 'hcrs', 'icrs', etc.,
+            or an actual coordinate system instance
 
         Returns
         -------

src/dysh/spectra/tests/test_doppler_frames.py

@@ -42,7 +42,7 @@ def test_compare_with_GBTIDL(self, data_dir):
             "HEL": "hcrs",
             "LSD": "lsrd",
             "LSR": "lsrk",
-            "TOPO": "topo",
+            "TOPO": sp.observer,  # must be an actual coordinate for topographical
         }
         maxHzdiff = {"BAR": 2.0, "GAL": 13000.0, "GEO": 2.0, "HEL": 2.0, "LSD": 30.0, "LSR": 2.0, "TOPO": 2.0}
         for k, v in framedict.items():

src/dysh/spectra/tests/test_spectrum.py

@@ -336,7 +336,7 @@ def test_smooth(self):
         for frame in astropy_frame_dict.keys():
             try:
                 s.set_frame(frame)
-            except KeyError:
+            except Exception:
                 print(f"set_frame fails for: {frame}")
                 continue
             print(f"frame set to: {frame}")
@@ -539,6 +539,59 @@ def test_average_spectra(self):
         compare_spectrum(ps0_org, self.ps0, ignore_history=True, ignore_comments=True)
         compare_spectrum(ps1_org, self.ps1, ignore_history=True, ignore_comments=True)
 
+    def test_spectrum_with_frame(self):
+        """Regression test for issue #401 to ensure Spectrum.with_frame functions as advertised.
+        https://github.com/GreenBankObservatory/dysh/issues/401
+        """
+        spec = Spectrum.fake_spectrum()
+        # Ensure that repeated changes of frame to the same frame do note
+        # change after the first transform
+        s1 = spec.with_frame("lsrk")
+        s2 = s1.with_frame("lsrk")
+        assert all(s1.spectral_axis == s2.spectral_axis)
+        assert s2.velocity_frame == "lsrk"
+        assert s2.meta["VELDEF"][4:] == "-LSR"
+
+        # Test that topographic results in an Exception because
+        # users must provide an ITRS coordinate instance in that case.
+        # First change the spectrum frame to somthing else because
+        # if it is already topo/itrs, then the error is circumvented
+        spec2 = Spectrum.fake_spectrum()
+        spec2.set_frame("gcrs")
+        with pytest.raises(ValueError):
+            spec2.set_frame("topo")
+
+        # Setting a new frame to the old frame does NOT result in an
+        # identical observer attribute on the resultant SpectralAxis.
+        # See https://github.com/astropy/astropy/issues/17506
+        # This test ensures that the difference remains small.
+
+        location_diff = np.sqrt(
+            (s1.observer.x - s2.observer.x) ** 2
+            + (s1.observer.y - s2.observer.y) ** 2
+            + (s1.observer.z - s2.observer.z) ** 2
+        )
+        velocity_diff = np.sqrt(
+            (s1.observer.v_x - s2.observer.v_x) ** 2
+            + (s1.observer.v_y - s2.observer.v_y) ** 2
+            + (s1.observer.v_z - s2.observer.v_z) ** 2
+        )
+        assert location_diff < 1.0e-5 * u.m
+        assert velocity_diff < 2e-8 * u.km / u.s
+
+    def test_velocity_axis_to(self):
+        """Regression test for issue 372 https://github.com/GreenBankObservatory/dysh/issues/372
+        Calling velocity_axis_to should not change the object spectral axis.
+        """
+        spec = Spectrum.fake_spectrum()
+        id1 = id(spec.spectral_axis)
+        sa = spec.velocity_axis_to(toframe="lsrk")
+        assert id(spec.spectral_axis) == id1
+        assert id(sa) != id1
+        # converting to identical frame should not change the values in the spectral axis
+        sa = spec.velocity_axis_to(toframe="topo")
+        assert all(sa == spec.spectral_axis)
+
     def test_baseline(self):
         """Test for comparing GBTIDL baseline to Dysh baselines"""