184 change the doppler convention

src/dysh/coordinates/__init__.py

@@ -2,3 +2,4 @@
 
 __all__ = ["core"]
 from .core import *
+from .core import KMS

src/dysh/coordinates/core.py

@@ -16,6 +16,7 @@
 _PMZERORAD = 0.0 * u.rad / u.s
 _VELZERO = 0.0 * u.km / u.s
 _MPS = u.m / u.s
+KMS = u.km / u.s
 
 # Velocity frame conventions, partially stolen from pyspeckit.
 # See also Section6.2.5.2 of the GBT observer's guide https://www.gb.nrao.edu/scienceDocs/GBTog.pdf
@@ -95,13 +96,19 @@
 # Dictionary to convert from FITS velocity convention to specutils string.
 # At GBT, VELO was written by sdfits filler for some unknown amount of
 # time instead of RELA, so allow for it here
-vconv_dict = {
+velocity_convention_dict = {
     "OPTI": "optical",
     "RADI": "radio",
     "RELA": "relativistic",
     "VELO": "relativistic",
 }
 
+reverse_velocity_convention_dict = {"optical": "OPTI", "radio": "RADI", "relativistic": "VELO"}
+
+
+def replace_convention(veldef, doppler_convention):
+    return reverse_velocity_convention_dict[doppler_convention] + veldef[4:]
+
 
 # This gives the wrong answer for GBT which always writes data as topocentric
 # regardless of what VELDEF says.
@@ -136,7 +143,7 @@ def decode_veldef(veldef):
         raise ValueError(f"VELDEF string {veldef} must be no more than 8 characters.")
     vconv = veldef[:4]
     try:
-        velocity_convention = vconv_dict[vconv]
+        velocity_convention = velocity_convention_dict[vconv]
     except KeyError:
         raise KeyError(f"Velocity convention {vconv} not recognized.")
 

src/dysh/spectra/spectrum.py

@@ -1,4 +1,5 @@
 import warnings
+from copy import deepcopy
 
 import astropy.units as u
 import numpy as np
@@ -11,16 +12,16 @@
 from astropy.wcs import WCS, FITSFixedWarning
 from specutils import Spectrum1D
 
-from ..coordinates import (
+from ..coordinates import (  # is_topocentric,; topocentric_velocity_to_frame,
+    KMS,
     Observatory,
     astropy_frame_dict,
     change_ctype,
     decode_veldef,
     get_velocity_in_frame,
-    is_topocentric,
     make_target,
+    replace_convention,
     sanitize_skycoord,
-    topocentric_velocity_to_frame,
     veldef_to_convention,
 )
 from ..plot import specplot as sp
@@ -49,6 +50,8 @@ def __init__(self, *args, **kwargs):
             self._velocity_frame = self._target.frame.name
         else:
             self._velocity_frame = None
+        # @TODO - have _observer_location attribute instead
+        # and observer property returns getITRS(observer_location,obstime)
         self._observer = kwargs.pop("observer", None)
         Spectrum1D.__init__(self, *args, **kwargs)
         self._spectral_axis._target = self._target
@@ -57,9 +60,11 @@ def __init__(self, *args, **kwargs):
             self._obstime = Time(self.meta["DATE-OBS"])
         else:
             self._obstime = None
-        if "CTYPE1" in self.meta:
-            self._target.topocentric = is_topocentric(self.meta["CTYPE1"])
-            self.topocentric = is_topocentric(self.meta["CTYPE1"])
+        # if "CTYPE1" in self.meta:
+        #    # may not need these attributes anymore.
+        #    if self._target is not None:
+        #        self._target.topocentric = is_topocentric(self.meta["CTYPE1"])
+        #    self.topocentric = is_topocentric(self.meta["CTYPE1"])
 
         # if mask is not set via the flux input (NaNs in flux or flux.mask),
         # then set the mask to all False == good
@@ -219,7 +224,7 @@ def equivalencies(self):
         else:
             rfq = None
         if rfq is not None:
-            equiv.extend(get_spectral_equivalency(rfq, self._velocity_convention))
+            equiv.extend(get_spectral_equivalency(rfq, self.velocity_convention))
         return equiv
 
     @property
@@ -234,43 +239,125 @@ def observer(self):
     def velocity_frame(self):
         return self._velocity_frame
 
-    @property
-    def velocity_convention(self):
-        return self._spectral_axis.doppler_convention
+    # This is already in specutils.OneDSpectrumMixin
+    # @property
+    # def velocity_convention(self):
+    #    return self._spectral_axis.doppler_convention
 
     @property
     def doppler_convention(self):
         return self.velocity_convention
 
-    def velocity_axis_to(self, unit=u.km / u.s, toframe=None, toconvention=None):
-        if toframe is not None:
-            self.shift_to_frame(toframe)
-        if toconvention is not None:
-            return self._spectral_axis.to(unit=unit, doppler_convention=toconvention)
-        else:
-            return self.velocity.to(unit)
+    def axis_velocity(self, unit=KMS):
+        """Get the spectral axis in velocity units.
+        *Note*: This is not the same as `Spectrum.velocity`, which includes the source radial velocity.
+        """
+        return self._spectral_axis.to(unit)
+
+    # not needed
+    # def velocity_axis_to(self, unit=KMS, toframe=None, doppler_convention=None):
+    #    if toframe is not None:
+    #        self.set_frame(toframe)
+    #    if doppler_convention is not None:
+    #        return self._spectral_axis.to(unit=unit, doppler_convention=doppler_convention).to(unit)
+    #    else:
+    #        return self.velocity.to(unit)
 
     def get_velocity_shift_to(self, toframe):
         if self._target is None:
             raise Exception("Can't calculate velocity because Spectrum.target is None")
         return get_velocity_in_frame(self._target, toframe, self._observer, self._obstime)
 
-    def shift_to_frame(self, toframe):
-        # v = self.get_velocity_shift_to(toframe)
-        # self._spectral_axis = self._spectral_axis.with_radial_velocity_shift(target_shift=v)
-        actualframe = astropy_frame_dict.get(toframe, toframe)
+    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.
+        """
+
+        if "topo" in toframe:
+            actualframe = self.observer  # ???
+        else:
+            actualframe = astropy_frame_dict.get(toframe, toframe)
         # print(f"actual frame is {actualframe}")
         self._spectral_axis = self._spectral_axis.with_observer_stationary_relative_to(actualframe)
         self._meta["CTYPE1"] = change_ctype(self._meta["CTYPE1"], toframe)
         # @todo shouldn't have two variables for the same thing.
         # Still needed?
-        self.topocentric = self._target.topocentric = "topo" in toframe
+        # self.topocentric = self._target.topocentric = "topo" in toframe
+
+    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.
+
+        Returns
+        -------
+        spectrum : `~dysh.spectra.Spectrum`
+            A new Spectrum object
+        """
+
+        s = deepcopy(self)
+        s.set_frame(toframe)
+        return s
+
+    def set_convention(self, doppler_convention):
+        """Set the velocity convention of this Spectrum.  The spectral axis of this Spectrum will be replaced
+        with a new spectral axis with the input velocity convention.  The header 'VELDEF' value will
+        be changed accordingly.
+
+        To make a copy of this Spectrum with a new velocity convention instead, use `with_velocity_convention`.
 
-    # Not sure we ever actually want to do this.
-    # The convention is a "display" parameter.
-    def _change_convention(self, toconvention):
-        new_sp_axis = self.spectral_axis.replicate(doppler_convention=toconvention)
+        Parameters
+        ----------
+        doppler_convention - str
+            The velocity convention, one of 'radio', 'optical', 'relativistic'
+
+        """
+        # replicate() gives the same asnwer as
+        # self._spectral_axis.to(unit=self._spectral_axis.unit, doppler_convention=doppler_convention)
+        new_sp_axis = self.spectral_axis.replicate(doppler_convention=doppler_convention)
         self._spectral_axis = new_sp_axis
+        self.meta["VELDEF"] = replace_convention(self.meta["VELDEF"], doppler_convention)
+
+    def with_velocity_convention(self, doppler_convention):
+        """Returns a copy of this Spectrum with the input velocity convention.  The header 'VELDEF' value will
+        be changed accordingly.
+
+        Parameters
+        ----------
+        doppler_convention - str
+            The velocity convention, one of 'radio', 'optical', 'relativistic'
+
+        Returns
+        -------
+        spectrum : `~dysh.spectra.Spectrum`
+            A new Spectrum object
+        """
+        if False:
+            # this doesn't work.
+            # for some reason, the axis velocity
+            # still contains the difference between TOPO and observed frame
+            s = self.__class__(
+                flux=self.flux,
+                wcs=self.wcs,
+                meta=self.meta,
+                velocity_convention=doppler_convention,
+                target=self._target,
+                observer=self._spectral_axis.observer,
+            )
+            s.meta["VELDEF"] = replace_convention(self.meta["VELDEF"], doppler_convention)
+        s = deepcopy(self)
+        s.set_convention(doppler_convention)
+        return s
 
     def savefig(self, file, **kwargs):
         """Save the plot"""
@@ -303,8 +390,36 @@ def _write_table(self, fileobj, format, **kwargs):
         # f=kwargs.pop("format")
         t.write(fileobj, format=format, **kwargs)
 
+    def _copy(self, **kwargs):
+        """
+        Perform deep copy operations on each attribute of the ``Spectrum``
+        object.
+        This overrides the ``specutils.Spectrum1D`` method so that
+        target and observer attributes get copied.
+        """
+        alt_kwargs = dict(
+            flux=deepcopy(self.flux),
+            spectral_axis=deepcopy(self.spectral_axis),
+            uncertainty=deepcopy(self.uncertainty),
+            wcs=deepcopy(self.wcs),
+            mask=deepcopy(self.mask),
+            meta=deepcopy(self.meta),
+            unit=deepcopy(self.unit),
+            velocity_convention=deepcopy(self.velocity_convention),
+            rest_value=deepcopy(self.rest_value),
+            target=deepcopy(self.target),
+            observer=deepcopy(self.observer),
+        )
+
+        alt_kwargs.update(kwargs)
+
+        s = self.__class__(**alt_kwargs)
+        # s.topocentric = is_topocentric(meta["CTYPE1"])  # GBT-specific to use CTYPE1 instead of VELDEF
+        return s
+
     @classmethod
-    def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, shift_topo=False):
+    def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None):
+        # , shift_topo=False):
         """Factory method to create a Spectrum object from a data and header.
 
         Parameters
@@ -359,7 +474,7 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, shift_t
             # reset warnings?
         else:
             wcs = None
-        is_topo = is_topocentric(meta["CTYPE1"])  # GBT-specific to use CTYPE1 instead of VELDEF
+        # is_topo = is_topocentric(meta["CTYPE1"])  # GBT-specific to use CTYPE1 instead of VELDEF
         target = make_target(meta)
         vc = veldef_to_convention(meta["VELDEF"])
         vf = astropy_frame_dict[decode_veldef(meta["VELDEF"])[1]]
@@ -397,9 +512,9 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, shift_t
         if observer_location is None:
             s.set_radial_velocity_to(target.radial_velocity)
         # I THINK THIS IS NO LONGER NEEDED
-        if shift_topo:  # and is_topo
-            vshift = topocentric_velocity_to_frame(target, vf, observer=Observatory["GBT"], obstime=obstime)
-
+        # if shift_topo:  # and is_topo
+        #    vshift = topocentric_velocity_to_frame(target, vf, observer=Observatory["GBT"], obstime=obstime)
+        #
         return s
 
 

src/dysh/spectra/tests/test_doppler_conventions.py

@@ -0,0 +1,41 @@
+from copy import deepcopy
+from gzip import GzipFile
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+from dysh.fits import gbtfitsload
+
+
+class TestDopplerConvention:
+    def read_ascii(self, filename, unzip=True):
+        if unzip:
+            g = GzipFile(filename)
+        else:
+            g = filename
+        values = np.loadtxt(g, skiprows=3, unpack=True)
+        velokms = values[:][0] * u.km / u.s
+        flux = values[:][1] * u.ct
+        return velokms, flux
+
+    def compare_gbtidl(self, spectrum, filename, doppler_convention, maxdiff):
+        velokms, flux = self.read_ascii(filename)
+        s = spectrum.with_velocity_convention(doppler_convention)
+        x = np.mean(velokms - s.axis_velocity())
+        assert np.abs(x.value) < maxdiff
+
+    def test_compare_with_GBTIDL(self, data_dir):
+        """ """
+        # get filenames
+        sdf_file = f"{data_dir}/TGBT21A_501_11/TGBT21A_501_11.raw.vegas.fits"
+
+        sdf = gbtfitsload.GBTFITSLoad(sdf_file)
+        sp = sdf.getps(152)[0].calibrated(0)
+        sp.set_frame("hcrs")  # data were tracked in heliocentric. change in place.
+
+        conventions = {"OPTI-HEL": "optical", "RADI-HEL": "radio", "TRUE-HEL": "relativistic"}
+        maxdiff = 1.0
+        for k, v in conventions.items():
+            gbtidl_file = f"{data_dir}/gbtidl_spectra/onoff-L_getps_152_{k}.ascii.gz"
+            self.compare_gbtidl(filename=gbtidl_file, spectrum=sp, doppler_convention=v, maxdiff=maxdiff)

src/dysh/spectra/tests/test_doppler_frames.py

@@ -1,13 +1,13 @@
-from copy import deepcopy
-
 import astropy.units as u
 import numpy as np
 import pytest
 
 from dysh.fits import gbtfitsload
 
 
-class TestDoppler:
+# @TODO consolidate methods with test_doppler_conventions,
+# possibly in util module
+class TestVelocityFrames:
     def read_ascii(self, filename):
         values = np.loadtxt(filename, skiprows=3, unpack=True)
         # print(np.shape(values))
@@ -17,8 +17,8 @@ def read_ascii(self, filename):
 
     def compare_gbtidl(self, spectrum, filename, frame, maxdiff):
         freqGHz, flux = self.read_ascii(filename)
-        s = deepcopy(spectrum)
-        s.shift_to_frame(frame)
+        # use version of method that makes a copy, so it gets some test coverage.
+        s = spectrum.with_frame(frame)
         x = np.mean(freqGHz - s.spectral_axis)
         y = np.nanmean(flux - s.flux)
         # print(x, x.to("Hz"))
@@ -42,15 +42,9 @@ def test_compare_with_GBTIDL(self, data_dir):
             "HEL": "hcrs",
             "LSD": "lsrd",
             "LSR": "lsrk",
-        }  #'TOPO':'topo'
-        maxHzdiff = {
-            "BAR": 2.0,
-            "GAL": 13000.0,
-            "GEO": 2.0,
-            "HEL": 2.0,
-            "LSD": 30.0,
-            "LSR": 2.0,
-        }  #'TOPO':'topo'
+            "TOPO": "topo",
+        }
+        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():
             gbtidl_file = f"{data_dir}/gbtidl_spectra/onoff-L_gettp_156_intnum_0_{k}.ascii"
             self.compare_gbtidl(filename=gbtidl_file, spectrum=sp, frame=v, maxdiff=maxHzdiff[k])