Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Spectrum1D translator update for NDCube version #36

Merged
merged 5 commits into from
Aug 27, 2021
Merged

Spectrum1D translator update for NDCube version #36

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
010b9ab
Update Spectrum1D translator to work with NDCube-based 3D cubes
rosteen Jul 8, 2021
2940e3d
Handling 1D dummy GWCS case
rosteen Jul 9, 2021
1002ca8
Working on better 3D Spectrum1D tests
rosteen Jul 19, 2021
cdd844d
Fix problems with uncertainty/mask axis swapping and collapsed coords
rosteen Jul 21, 2021
7a01dbf
Working on new getting new unit test to pass
rosteen Jul 21, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 31 additions & 10 deletions glue_astronomy/translators/spectrum1d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -21,26 +22,44 @@
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).
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In addition to the above comments it would be good to clarify here whether you mean for the WCS or Numpy array. The data array and WCS will have axes in opposite order - spectral cubes read with plain Astropy FITS and WCS or with SpectralCube are often stored in spectral,dec,ra for the Numpy data array and ra,dec,spectral for the WCS.

Glue can handle any axis order as long as x_att and y_att are set correctly but it would indeed be good to make sure SpectralCube and Spectrum1D return a similar glue data object.

# 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:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about cases where the data has 2 dimensions? (for example a 2D spectrum)

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:
# 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))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't really need this anymore, Data coords can be a GWCS as far as I know.

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(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(obj.mask, -1, 0)
else:
data['mask'] = obj.mask

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.

Expand Down Expand Up @@ -69,14 +88,16 @@ def to_object(self, data_or_subset, attribute=None, statistic='mean'):
# Find non-spectral axes
axes = tuple(i for i in range(data.ndim) if i != spec_axis)

kwargs = {'wcs': data.coords.sub([WCSSUB_SPECTRAL])}
if statistic is not None:
kwargs = {'wcs': data.coords.sub([WCSSUB_SPECTRAL])}
else:
kwargs = {'wcs': data.coords}

elif isinstance(data.coords, SpectralCoordinates):

kwargs = {'spectral_axis': data.coords.spectral_axis}

else:

raise TypeError('data.coords should be an instance of WCS '
'or SpectralCoordinates')

Expand Down Expand Up @@ -112,7 +133,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)
Expand Down
62 changes: 43 additions & 19 deletions glue_astronomy/translators/tests/test_spectrum1d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,23 +123,36 @@ 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':
# 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', 'FREQ']
wcs.wcs.set()
kwargs = {'wcs': wcs}
flux = np.ones((4, 4, 5))*u.Unit('Jy')
uncertainty = VarianceUncertainty(np.square(flux*0.1))
mask = np.zeros((4, 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()

Expand All @@ -150,27 +163,38 @@ 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'

# Check round-tripping via single attribute reference
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)