From 677650501c3748edc470626af54f2e60ac4b1436 Mon Sep 17 00:00:00 2001 From: "busko@stsci.edu" Date: Wed, 24 Jul 2019 11:16:08 -0400 Subject: [PATCH 1/7] Add support for redshift in conversions to velocity. --- specutils/spectra/spectrum_mixin.py | 7 +++-- specutils/tests/test_spectrum1d.py | 42 +++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py index 0b3da9253..8d1216507 100644 --- a/specutils/spectra/spectrum_mixin.py +++ b/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,9 @@ def velocity(self): new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv) + # Redshift-correct velocities using approximate low-z correction. + new_data += Quantity((self._redshift * C), new_data.unit) + return new_data def with_spectral_unit(self, unit, velocity_convention=None, diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py index 510dfbefa..213e603fa 100644 --- a/specutils/tests/test_spectrum1d.py +++ b/specutils/tests/test_spectrum1d.py @@ -106,6 +106,48 @@ def test_spectral_axis_conversions(): new_spec = spec.with_spectral_unit(u.GHz) +def test_redshift(): + spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, + spectral_axis=np.array([4000, 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, 0.5) + assert spec.velocity[1].value == pytest.approx(99930.8, 0.5) + + spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, + spectral_axis=np.array([4000, 8000]) * u.AA, + velocity_convention='optical', + rest_value=6000 * u.AA, + redshift = 0.1) + + assert spec.velocity.unit == u.Unit('km/s') + assert spec.velocity[0].value == pytest.approx(-69930.8, 0.5) + assert spec.velocity[1].value == pytest.approx(129930.8, 0.5) + + #------------------------- + + spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, + spectral_axis=np.array([10.5, 11.5]) * u.GHz, + velocity_convention='radio', + rest_value=11. * u.GHz) + + assert spec.velocity.unit == u.Unit('km/s') + assert spec.velocity[0].value == pytest.approx(13626., 1.) + assert spec.velocity[1].value == pytest.approx(-13626., 1.) + + spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, + spectral_axis=np.array([10.5, 11.5]) * u.GHz, + velocity_convention='radio', + rest_value=11. * u.GHz, + redshift = 0.1) + + assert spec.velocity.unit == u.Unit('km/s') + assert spec.velocity[0].value == pytest.approx(43606., 1.) + assert spec.velocity[1].value == pytest.approx(16352., 1.) + + 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, From 2bf680ac45372f0a0b8d5e6ef776610e17ffa10d Mon Sep 17 00:00:00 2001 From: "busko@stsci.edu" Date: Wed, 24 Jul 2019 11:17:02 -0400 Subject: [PATCH 2/7] Add support for redshift in conversions to velocity. --- setup.cfg | 2 +- specutils/spectra/spectrum1d.py | 9 ++++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/setup.cfg b/setup.cfg index 86950319e..27a2ec80a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,7 +21,7 @@ remote_data_strict = True asdf_schema_root = specutils/io/asdf/schemas # The remote data tests will run by default. Passing --remote-data=none on the # command line will override this setting. -addopts = --remote-data=any --doctest-rst +#addopts = --remote-data=any --doctest-rst [ah_bootstrap] auto_use = True diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py index 08d88312b..0ac47126f 100644 --- a/specutils/spectra/spectrum1d.py +++ b/specutils/spectra/spectrum1d.py @@ -32,6 +32,9 @@ 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 : float + The redshift correction to be applied to the spectral axis for use + with velocity conversions. 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 +44,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=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 +111,10 @@ 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.") + self._redshift = redshift + if redshift is None: + self._redshift = 0. + super(Spectrum1D, self).__init__( data=flux.value if isinstance(flux, u.Quantity) else flux, wcs=wcs, **kwargs) From 22928f9bf1d4e501476a3eda2cc1220f846f079d Mon Sep 17 00:00:00 2001 From: "busko@stsci.edu" Date: Wed, 24 Jul 2019 11:19:09 -0400 Subject: [PATCH 3/7] Revert change in config file. --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 27a2ec80a..86950319e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,7 +21,7 @@ remote_data_strict = True asdf_schema_root = specutils/io/asdf/schemas # The remote data tests will run by default. Passing --remote-data=none on the # command line will override this setting. -#addopts = --remote-data=any --doctest-rst +addopts = --remote-data=any --doctest-rst [ah_bootstrap] auto_use = True From 0b6d17c25aed2b3a7c716161c4bf2c74f7fd793c Mon Sep 17 00:00:00 2001 From: "busko@stsci.edu" Date: Wed, 24 Jul 2019 11:58:08 -0400 Subject: [PATCH 4/7] Fine tune tests. --- specutils/tests/test_spectrum1d.py | 40 ++++++++++++++++-------------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py index 213e603fa..436cb38ef 100644 --- a/specutils/tests/test_spectrum1d.py +++ b/specutils/tests/test_spectrum1d.py @@ -107,45 +107,49 @@ def test_spectral_axis_conversions(): def test_redshift(): - spec = Spectrum1D(flux=np.array([26.0, 44.5]) * u.Jy, - spectral_axis=np.array([4000, 8000]) * u.AA, + 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, 0.5) - assert spec.velocity[1].value == pytest.approx(99930.8, 0.5) + 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, 44.5]) * u.Jy, - spectral_axis=np.array([4000, 8000]) * u.AA, + 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 = 0.1) assert spec.velocity.unit == u.Unit('km/s') - assert spec.velocity[0].value == pytest.approx(-69930.8, 0.5) - assert spec.velocity[1].value == pytest.approx(129930.8, 0.5) + 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, 44.5]) * u.Jy, - spectral_axis=np.array([10.5, 11.5]) * u.GHz, + 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. * u.GHz) + rest_value=11.0 * u.GHz) assert spec.velocity.unit == u.Unit('km/s') - assert spec.velocity[0].value == pytest.approx(13626., 1.) - assert spec.velocity[1].value == pytest.approx(-13626., 1.) + 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, 44.5]) * u.Jy, - spectral_axis=np.array([10.5, 11.5]) * u.GHz, + 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. * u.GHz, + rest_value=11.0 * u.GHz, redshift = 0.1) assert spec.velocity.unit == u.Unit('km/s') - assert spec.velocity[0].value == pytest.approx(43606., 1.) - assert spec.velocity[1].value == pytest.approx(16352., 1.) + 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.) def test_flux_unit_conversion(): From b1fb951a2535ca366431251314d2f1ba6dac960b Mon Sep 17 00:00:00 2001 From: busko Date: Wed, 7 Aug 2019 11:36:45 -0400 Subject: [PATCH 5/7] Add support for radial velocity. --- specutils/spectra/spectrum1d.py | 16 +++++++++------- specutils/spectra/spectrum_mixin.py | 5 ++++- specutils/tests/test_spectrum1d.py | 28 ++++++++++++++++++++++++++-- 3 files changed, 39 insertions(+), 10 deletions(-) diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py index 0ac47126f..b408b9f7f 100644 --- a/specutils/spectra/spectrum1d.py +++ b/specutils/spectra/spectrum1d.py @@ -32,9 +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 : float - The redshift correction to be applied to the spectral axis for use - with velocity conversions. + redshift : 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 @@ -44,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, redshift=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): @@ -111,9 +113,9 @@ 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.") - self._redshift = redshift - if redshift is None: - self._redshift = 0. + 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, diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py index 8d1216507..ea46fc020 100644 --- a/specutils/spectra/spectrum_mixin.py +++ b/specutils/spectra/spectrum_mixin.py @@ -186,7 +186,10 @@ def velocity(self): new_data = self.spectral_axis.to(u.km/u.s, equivalencies=equiv) # Redshift-correct velocities using approximate low-z correction. - new_data += Quantity((self._redshift * C), new_data.unit) + 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 diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py index 436cb38ef..b945964a3 100644 --- a/specutils/tests/test_spectrum1d.py +++ b/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 @@ -121,7 +122,7 @@ def test_redshift(): spectral_axis=np.array([4000, 6000, 8000]) * u.AA, velocity_convention='optical', rest_value=6000 * u.AA, - redshift = 0.1) + redshift_rv= 0.1) assert spec.velocity.unit == u.Unit('km/s') assert spec.velocity[0].value == pytest.approx(-69951.3, abs=0.5) @@ -144,13 +145,36 @@ def test_redshift(): spectral_axis=np.array([10.5, 11.0, 11.5]) * u.GHz, velocity_convention='radio', rest_value=11.0 * u.GHz, - redshift = 0.1) + 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 From ebff93442cf26d287302f4747a1663d152306d24 Mon Sep 17 00:00:00 2001 From: Ivo Busko Date: Fri, 9 Aug 2019 11:03:25 -0400 Subject: [PATCH 6/7] Fix comment and doc string --- specutils/spectra/spectrum1d.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py index b408b9f7f..d00ef7385 100644 --- a/specutils/spectra/spectrum1d.py +++ b/specutils/spectra/spectrum1d.py @@ -32,7 +32,7 @@ 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 : float or `~astropy.units.Quantity` + 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, @@ -113,6 +113,8 @@ 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. From 44751c19b3b2873e5671f5af3d3b90e62e496cbc Mon Sep 17 00:00:00 2001 From: Ivo Busko Date: Fri, 9 Aug 2019 11:38:05 -0400 Subject: [PATCH 7/7] Attempt to fix whitespace --- specutils/spectra/spectrum1d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/specutils/spectra/spectrum1d.py b/specutils/spectra/spectrum1d.py index d00ef7385..29e995f02 100644 --- a/specutils/spectra/spectrum1d.py +++ b/specutils/spectra/spectrum1d.py @@ -118,7 +118,7 @@ def __init__(self, flux=None, spectral_axis=None, wcs=None, 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)