From 010b9abbbde767e278d7b897ee1ec548bd62d59d Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Thu, 8 Jul 2021 17:46:26 -0400 Subject: [PATCH 1/5] Update Spectrum1D translator to work with NDCube-based 3D cubes --- glue_astronomy/translators/spectrum1d.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/glue_astronomy/translators/spectrum1d.py b/glue_astronomy/translators/spectrum1d.py index 42c22a9..9e19e76 100644 --- a/glue_astronomy/translators/spectrum1d.py +++ b/glue_astronomy/translators/spectrum1d.py @@ -21,26 +21,38 @@ class Specutils1DHandler: def to_data(self, obj): - coords = SpectralCoordinates(obj.spectral_axis) - data = Data(coords=coords) - data['flux'] = obj.flux - data.get_component('flux').units = str(obj.unit) + + # Glue expects spectral axis first for cubes (opposite of specutils). + # Swap the spectral axis to first here. to_object doesn't need this because + # Spectrum1D does it automatically on initialization. + if len(obj.flux.shape) == 3: + data = Data(coords=obj.wcs.swapaxes(-1, 0)) + data['flux'] = np.swapaxes(obj.flux, -1, 0) + data.get_component('flux').units = str(obj.unit) + else: + data = Data(coords=obj.wcs) + data['flux'] = obj.flux + data.get_component('flux').units = str(obj.unit) # Include uncertainties if they exist if obj.uncertainty is not None: data['uncertainty'] = obj.uncertainty.quantity + if len(obj.flux.shape) == 3: + data['uncertainty'] = np.swapaxes(data['uncertainty'], -1, 0) data.get_component('uncertainty').units = str(obj.uncertainty.unit) data.meta.update({'uncertainty_type': obj.uncertainty.uncertainty_type}) # Include mask if it exists if obj.mask is not None: data['mask'] = obj.mask + if len(obj.flux.shape) == 3: + data['mask'] = np.swapaxes(data['mask'], -1, 0) data.meta.update(obj.meta) return data - def to_object(self, data_or_subset, attribute=None, statistic='mean'): + def to_object(self, data_or_subset, attribute=None, statistic=None): """ Convert a glue Data object to a Spectrum1D object. @@ -112,7 +124,7 @@ def parse_attributes(attributes): mask = ~mask # Collapse values and mask to profile - if data.ndim > 1: + if data.ndim > 1 and statistic is not None: # Get units and attach to value values = data.compute_statistic(statistic, attribute, axis=axes, subset_state=subset_state) From 2940e3d0dbe7bcf29f07f2150f3fb5eadcb1a0d0 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Fri, 9 Jul 2021 13:17:10 -0400 Subject: [PATCH 2/5] Handling 1D dummy GWCS case Fix to use SpectralCoordinate --- glue_astronomy/translators/spectrum1d.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/glue_astronomy/translators/spectrum1d.py b/glue_astronomy/translators/spectrum1d.py index 9e19e76..5b4fe93 100644 --- a/glue_astronomy/translators/spectrum1d.py +++ b/glue_astronomy/translators/spectrum1d.py @@ -7,6 +7,7 @@ from astropy import units as u from astropy.wcs import WCSSUB_SPECTRAL from astropy.nddata import StdDevUncertainty, InverseVariance, VarianceUncertainty +from gwcs import WCS as GWCS from glue_astronomy.spectral_coordinates import SpectralCoordinates @@ -30,7 +31,11 @@ def to_data(self, obj): data['flux'] = np.swapaxes(obj.flux, -1, 0) data.get_component('flux').units = str(obj.unit) else: - data = Data(coords=obj.wcs) + # Don't use the dummy GWCS created by Spectrum1D initialized with spectral_axis + if isinstance(obj.wcs, GWCS): + data = Data(coords=SpectralCoordinates(obj.spectral_axis)) + else: + data = Data(coords=obj.wcs) data['flux'] = obj.flux data.get_component('flux').units = str(obj.unit) @@ -88,7 +93,7 @@ def to_object(self, data_or_subset, attribute=None, statistic=None): kwargs = {'spectral_axis': data.coords.spectral_axis} else: - + print(type(data.coords)) raise TypeError('data.coords should be an instance of WCS ' 'or SpectralCoordinates') From 1002ca89ecbc469c006754fbb4dfb341414ce8f4 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Mon, 19 Jul 2021 16:08:16 -0400 Subject: [PATCH 3/5] Working on better 3D Spectrum1D tests Working on better 3D Spectrum1D tests --- glue_astronomy/translators/spectrum1d.py | 1 + .../translators/tests/test_spectrum1d.py | 44 ++++++++++++------- 2 files changed, 30 insertions(+), 15 deletions(-) diff --git a/glue_astronomy/translators/spectrum1d.py b/glue_astronomy/translators/spectrum1d.py index 5b4fe93..46818a8 100644 --- a/glue_astronomy/translators/spectrum1d.py +++ b/glue_astronomy/translators/spectrum1d.py @@ -87,6 +87,7 @@ def to_object(self, data_or_subset, attribute=None, statistic=None): axes = tuple(i for i in range(data.ndim) if i != spec_axis) kwargs = {'wcs': data.coords.sub([WCSSUB_SPECTRAL])} + #kwargs = {'wcs': data.coords} elif isinstance(data.coords, SpectralCoordinates): diff --git a/glue_astronomy/translators/tests/test_spectrum1d.py b/glue_astronomy/translators/tests/test_spectrum1d.py index 4ccfc61..162fe6a 100644 --- a/glue_astronomy/translators/tests/test_spectrum1d.py +++ b/glue_astronomy/translators/tests/test_spectrum1d.py @@ -72,6 +72,7 @@ def test_to_spectrum1d_from_3d_cube(): # Set up simple spectral WCS wcs = WCS(naxis=3) wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'VELO-LSR'] + #wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE-LOG'] wcs.wcs.set() data = Data(label='spectral-cube', coords=wcs) @@ -123,23 +124,35 @@ def test_to_spectrum1d_default_attribute(): 'keyword argument.') -@pytest.mark.parametrize('mode', ('wcs', 'lookup')) +@pytest.mark.parametrize('mode', ('wcs1d', 'wcs3d', 'lookup')) def test_from_spectrum1d(mode): - if mode == 'wcs': - wcs = WCS(naxis=1) - wcs.wcs.ctype = ['FREQ'] + if mode == 'wcs3d': + # Set up simple spectral WCS + wcs = WCS(naxis=3) + wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE-LOG'] + wcs.wcs.crval = [205, 27, 3.6216e-07] + wcs.wcs.crpix = [1, 1, 1] + wcs.wcs.cdelt = [-1.4e-4, 1.4e-4, 8.34e-11] wcs.wcs.set() - kwargs = {'wcs': wcs} + flux = np.ones((3, 4, 5))*u.Unit('Jy') + uncertainty = VarianceUncertainty(np.square(flux*0.1)) + mask = np.zeros((3,4,5)) + kwargs = {'wcs': wcs, 'uncertainty': uncertainty, 'mask': mask} else: - - kwargs = {'spectral_axis': [1, 2, 3, 4] * u.Hz} - - spec = Spectrum1D([2, 3, 4, 5] * u.Jy, - uncertainty=VarianceUncertainty( - [0.1, 0.1, 0.1, 0.1] * u.Jy**2), - mask=[False, False, False, False], - **kwargs) + flux = [2, 3, 4, 5] * u.Jy + uncertainty = VarianceUncertainty([0.1, 0.1, 0.1, 0.1] * u.Jy**2) + mask=[False, False, False, False] + if mode == 'wcs1d': + wcs = WCS(naxis=1) + wcs.wcs.ctype = ['FREQ'] + wcs.wcs.set() + kwargs = {'wcs': wcs, 'uncertainty': uncertainty, 'mask': mask} + else: + kwargs = {'spectral_axis': [1, 2, 3, 4] * u.Hz, + 'uncertainty': uncertainty, 'mask': mask} + + spec = Spectrum1D(flux, **kwargs) data_collection = DataCollection() @@ -150,13 +163,13 @@ def test_from_spectrum1d(mode): assert isinstance(data, Data) assert len(data.main_components) == 3 assert data.main_components[0].label == 'flux' - assert_allclose(data['flux'], [2, 3, 4, 5]) + assert_allclose(data['flux'], flux.value) component = data.get_component('flux') assert component.units == 'Jy' # Check uncertainty parsing within glue data object assert data.main_components[1].label == 'uncertainty' - assert_allclose(data['uncertainty'], [0.1, 0.1, 0.1, 0.1]) + assert_allclose(data['uncertainty'], uncertainty.array) component = data.get_component('uncertainty') assert component.units == 'Jy2' @@ -174,3 +187,4 @@ def test_from_spectrum1d(mode): assert_quantity_allclose(spec_new.flux, [2, 3, 4, 5] * u.Jy) assert spec_new.uncertainty is not None assert_quantity_allclose(spec_new.uncertainty.quantity, [0.1, 0.1, 0.1, 0.1] * u.Jy**2) + From cdd844d9cce08b30e4318d2e5110ab7bb2b44ba9 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 21 Jul 2021 15:50:08 -0400 Subject: [PATCH 4/5] Fix problems with uncertainty/mask axis swapping and collapsed coords --- glue_astronomy/translators/spectrum1d.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/glue_astronomy/translators/spectrum1d.py b/glue_astronomy/translators/spectrum1d.py index 46818a8..0245de2 100644 --- a/glue_astronomy/translators/spectrum1d.py +++ b/glue_astronomy/translators/spectrum1d.py @@ -41,17 +41,19 @@ def to_data(self, obj): # Include uncertainties if they exist if obj.uncertainty is not None: - data['uncertainty'] = obj.uncertainty.quantity if len(obj.flux.shape) == 3: - data['uncertainty'] = np.swapaxes(data['uncertainty'], -1, 0) + data['uncertainty'] = np.swapaxes(obj.uncertainty.quantity, -1, 0) + else: + data['uncertainty'] = obj.uncertainty.quantity data.get_component('uncertainty').units = str(obj.uncertainty.unit) data.meta.update({'uncertainty_type': obj.uncertainty.uncertainty_type}) # Include mask if it exists if obj.mask is not None: - data['mask'] = obj.mask if len(obj.flux.shape) == 3: - data['mask'] = np.swapaxes(data['mask'], -1, 0) + data['mask'] = np.swapaxes(obj.mask, -1, 0) + else: + data['mask'] = obj.mask data.meta.update(obj.meta) @@ -86,8 +88,10 @@ def to_object(self, data_or_subset, attribute=None, statistic=None): # Find non-spectral axes axes = tuple(i for i in range(data.ndim) if i != spec_axis) - kwargs = {'wcs': data.coords.sub([WCSSUB_SPECTRAL])} - #kwargs = {'wcs': data.coords} + if statistic is not None: + kwargs = {'wcs': data.coords.sub([WCSSUB_SPECTRAL])} + else: + kwargs = {'wcs': data.coords} elif isinstance(data.coords, SpectralCoordinates): From 7a01dbf0c552b2e5bc7eaea7dda51d83ef4b47c1 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 21 Jul 2021 16:58:07 -0400 Subject: [PATCH 5/5] Working on new getting new unit test to pass Finish setting up new Spectrum1D tests Fix codestyle errors Allow test to run on 2.0.dev --- glue_astronomy/translators/spectrum1d.py | 1 - .../translators/tests/test_spectrum1d.py | 38 ++++++++++++------- 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/glue_astronomy/translators/spectrum1d.py b/glue_astronomy/translators/spectrum1d.py index 0245de2..e2fa102 100644 --- a/glue_astronomy/translators/spectrum1d.py +++ b/glue_astronomy/translators/spectrum1d.py @@ -98,7 +98,6 @@ def to_object(self, data_or_subset, attribute=None, statistic=None): kwargs = {'spectral_axis': data.coords.spectral_axis} else: - print(type(data.coords)) raise TypeError('data.coords should be an instance of WCS ' 'or SpectralCoordinates') diff --git a/glue_astronomy/translators/tests/test_spectrum1d.py b/glue_astronomy/translators/tests/test_spectrum1d.py index 162fe6a..357f4f5 100644 --- a/glue_astronomy/translators/tests/test_spectrum1d.py +++ b/glue_astronomy/translators/tests/test_spectrum1d.py @@ -72,7 +72,6 @@ def test_to_spectrum1d_from_3d_cube(): # Set up simple spectral WCS wcs = WCS(naxis=3) wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'VELO-LSR'] - #wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE-LOG'] wcs.wcs.set() data = Data(label='spectral-cube', coords=wcs) @@ -128,21 +127,22 @@ def test_to_spectrum1d_default_attribute(): def test_from_spectrum1d(mode): if mode == 'wcs3d': - # Set up simple spectral WCS + # This test is intended to be run with the version of Spectrum1D based + # on NDCube 2.0 + pytest.importorskip("ndcube", minversion="1.99") + + # Set up simple spatial+spectral WCS wcs = WCS(naxis=3) - wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'WAVE-LOG'] - wcs.wcs.crval = [205, 27, 3.6216e-07] - wcs.wcs.crpix = [1, 1, 1] - wcs.wcs.cdelt = [-1.4e-4, 1.4e-4, 8.34e-11] + wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN', 'FREQ'] wcs.wcs.set() - flux = np.ones((3, 4, 5))*u.Unit('Jy') + flux = np.ones((4, 4, 5))*u.Unit('Jy') uncertainty = VarianceUncertainty(np.square(flux*0.1)) - mask = np.zeros((3,4,5)) + mask = np.zeros((4, 4, 5)) kwargs = {'wcs': wcs, 'uncertainty': uncertainty, 'mask': mask} else: flux = [2, 3, 4, 5] * u.Jy uncertainty = VarianceUncertainty([0.1, 0.1, 0.1, 0.1] * u.Jy**2) - mask=[False, False, False, False] + mask = [False, False, False, False] if mode == 'wcs1d': wcs = WCS(naxis=1) wcs.wcs.ctype = ['FREQ'] @@ -177,14 +177,24 @@ def test_from_spectrum1d(mode): spec_new = data.get_object(attribute='flux') assert isinstance(spec_new, Spectrum1D) assert_quantity_allclose(spec_new.spectral_axis, [1, 2, 3, 4] * u.Hz) - assert_quantity_allclose(spec_new.flux, [2, 3, 4, 5] * u.Jy) + if mode == 'wcs3d': + assert_quantity_allclose(spec_new.flux, np.ones((5, 4, 4))*u.Unit('Jy')) + else: + assert_quantity_allclose(spec_new.flux, [2, 3, 4, 5] * u.Jy) assert spec_new.uncertainty is None # Check complete round-tripping, including uncertainties spec_new = data.get_object() assert isinstance(spec_new, Spectrum1D) assert_quantity_allclose(spec_new.spectral_axis, [1, 2, 3, 4] * u.Hz) - assert_quantity_allclose(spec_new.flux, [2, 3, 4, 5] * u.Jy) - assert spec_new.uncertainty is not None - assert_quantity_allclose(spec_new.uncertainty.quantity, [0.1, 0.1, 0.1, 0.1] * u.Jy**2) - + if mode == 'wcs3d': + assert_quantity_allclose(spec_new.flux, np.ones((5, 4, 4))*u.Unit('Jy')) + assert spec_new.uncertainty is not None + print(spec_new.uncertainty) + print(uncertainty) + assert_quantity_allclose(spec_new.uncertainty.quantity, + np.ones((5, 4, 4))*0.01*u.Jy**2) + else: + assert_quantity_allclose(spec_new.flux, [2, 3, 4, 5] * u.Jy) + assert spec_new.uncertainty is not None + assert_quantity_allclose(spec_new.uncertainty.quantity, [0.1, 0.1, 0.1, 0.1] * u.Jy**2)