Merge 44751c1 into 7eaeb57

specutils/spectra/spectrum1d.py

@@ -32,6 +32,11 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         Any quantity supported by the standard spectral equivalencies
         (wavelength, energy, frequency, wave number). Describes the rest value
         of the spectral axis for use with velocity conversions.
+    redshift_rv : float or `~astropy.units.Quantity`
+        The redshift/RV correction to be applied to the spectral axis for use
+        with velocity conversions. If float, the value is interpreted as the
+        redshift z. If an instance of `~astropy.units.Quantity` is passed,
+        it is interpreted as a radial velocity with the appropriate velocity units.
     uncertainty : `~astropy.nddata.NDUncertainty`
         Contains uncertainty information along with propagation rules for
         spectrum arithmetic. Can take a unit, but if none is given, will use
@@ -41,7 +46,7 @@ class Spectrum1D(OneDSpectrumMixin, NDDataRef):
         around with the spectrum container object.
     """
     def __init__(self, flux=None, spectral_axis=None, wcs=None,
-                 velocity_convention=None, rest_value=None, *args, **kwargs):
+                 velocity_convention=None, rest_value=None, redshift_rv=None, *args, **kwargs):
         # If the flux (data) argument is a subclass of nddataref (as it would
         # be for internal arithmetic operations), avoid setup entirely.
         if isinstance(flux, NDDataRef):
@@ -108,6 +113,12 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None,
             elif not self._rest_value.unit.is_equivalent(u.AA) and not self._rest_value.unit.is_equivalent(u.Hz):
                 raise u.UnitsError("Rest value must be energy/wavelength/frequency equivalent.")
 
+        # redshift may be moved to somewhere else in the future,
+        # perhaps to the object that handles coordinates. 
+        self._redshift_rv = redshift_rv
+        if redshift_rv is None:
+            self._redshift_rv = 0.
+
         super(Spectrum1D, self).__init__(
             data=flux.value if isinstance(flux, u.Quantity) else flux,
             wcs=wcs, **kwargs)

specutils/spectra/spectrum_mixin.py

@@ -2,8 +2,8 @@
 from copy import deepcopy
 
 import numpy as np
-from astropy.wcs import WCSSUB_SPECTRAL
-from astropy.units import Unit
+from astropy.units import Quantity
+from astropy.constants import c as C
 from astropy import units as u
 import astropy.units.equivalencies as eq
 from astropy.utils.decorators import lazyproperty
@@ -185,6 +185,12 @@ def velocity(self):
 
         new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv)
 
+        # Redshift-correct velocities using approximate low-z correction.
+        if isinstance(self._redshift_rv, Quantity):
+            new_data += self._redshift_rv
+        else:
+            new_data += Quantity((self._redshift_rv * C), new_data.unit)
+
         return new_data
 
     def with_spectral_unit(self, unit, velocity_convention=None,

specutils/tests/test_spectrum1d.py

@@ -3,6 +3,7 @@
 import gwcs
 import numpy as np
 import pytest
+from astropy.units import Quantity
 from astropy.nddata import StdDevUncertainty
 
 from .conftest import remote_access
@@ -106,6 +107,75 @@ def test_spectral_axis_conversions():
     new_spec = spec.with_spectral_unit(u.GHz)
 
 
+def test_redshift():
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
+
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA,
+                      redshift_rv= 0.1)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-69951.3, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(29979.2, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(129910.1, abs=0.5)
+
+    #-------------------------
+
+    spec = Spectrum1D(flux=np.array([26.0, 30.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
+                      velocity_convention='radio',
+                      rest_value=11.0 * u.GHz)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(13626., abs=1.)
+    assert spec.velocity[1].value == pytest.approx(0., abs=1.)
+    assert spec.velocity[2].value == pytest.approx(-13626., abs=1.)
+
+    spec = Spectrum1D(flux=np.array([26.0, 30.0, 44.5]) * u.Jy,
+                      spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz,
+                      velocity_convention='radio',
+                      rest_value=11.0 * u.GHz,
+                      redshift_rv= 0.1)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(43606., abs=1.)
+    assert spec.velocity[1].value == pytest.approx(29979., abs=1.)
+    assert spec.velocity[2].value == pytest.approx(16352., abs=1.)
+
+    #------------------------- radial velocity mode
+
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA)
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-99930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(0.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(99930.8, abs=0.5)
+
+    spec = Spectrum1D(flux=np.array([26.0, 30., 44.5]) * u.Jy,
+                      spectral_axis=np.array([4000, 6000, 8000]) * u.AA,
+                      velocity_convention='optical',
+                      rest_value=6000 * u.AA,
+                      redshift_rv=Quantity(1000., 'km/s'))
+
+    assert spec.velocity.unit == u.Unit('km/s')
+    assert spec.velocity[0].value == pytest.approx(-98930.8, abs=0.5)
+    assert spec.velocity[1].value == pytest.approx(1000.0, abs=0.5)
+    assert spec.velocity[2].value == pytest.approx(100930.8, abs=0.5)
+
+
 def test_flux_unit_conversion():
     # By default the flux units should be set to Jy
     s = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy,