Update Mosviz to use Spectrum1D from SpectralCube

docs/mosviz/plugins.rst

@@ -69,4 +69,4 @@ The :guilabel:`Remove` button can be used to remove a slit once it has been appl
 In order to plot a slit onto the image viewer, we need WCS information from an image and slit position from a 2D spectrum.
 WCS information is taken from the `meta` attribute of the :class:`~astropy.nddata.CCDData` object representing the data in the
 image viewer. The slit position is calculated using the `S_REGION` header extension value, located in the `meta` attribute of
-the :class:`~spectral_cube.SpectralCube` object that is active in the 2D spectrum viewer.
+the :class:`~specutils.Spectrum1D` object that is active in the 2D spectrum viewer.

jdaviz/configs/mosviz/plugins/parsers.py

@@ -5,12 +5,12 @@
 from pathlib import Path
 
 import numpy as np
-from asdf.fits_embed import AsdfInFits
+from astropy import units as u
 from astropy.io import fits
+from astropy.io.registry import IORegistryError
 from astropy.wcs import WCS
 from glue.core.data import Data
 from glue.core.link_helpers import LinkSame
-from spectral_cube import SpectralCube
 from specutils import Spectrum1D, SpectrumList, SpectrumCollection
 
 from jdaviz.configs.imviz.plugins.parsers import get_image_data_iterator
@@ -110,14 +110,29 @@ def link_data_in_table(app, data_obj=None):
     # context manager, b) handling intra-row linkage of 1D and 2D spectra in a
     # loop, and c) handling inter-row linkage after that in one fell swoop.
     with app.data_collection.delay_link_manager_update():
-        for index in range(len(mos_data.get_component('1D Spectra').data)):
-            spec_1d = mos_data.get_component('1D Spectra').data[index]
-            spec_2d = mos_data.get_component('2D Spectra').data[index]
+
+        spectra_1d = mos_data.get_component('1D Spectra').data
+        spectra_2d = mos_data.get_component('2D Spectra').data
+
+        # Link each 1D spectrum with its corresponding 2D spectrum
+        for index in range(len(spectra_1d)):
+
+            spec_1d = spectra_1d[index]
+            spec_2d = spectra_2d[index]
 
             wc_spec_1d = app.session.data_collection[spec_1d].world_component_ids
             wc_spec_2d = app.session.data_collection[spec_2d].world_component_ids
 
-            wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_spec_2d[0]))
+            wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_spec_2d[1]))
+
+        # Link each 1D spectrum to all other 1D spectra
+        first_spec_1d = spectra_1d[0]
+        wc_first_spec_1d = app.session.data_collection[first_spec_1d].world_component_ids
+
+        for index in range(1, len(spectra_1d)):
+            spec_1d = spectra_1d[index]
+            wc_spec_1d = app.session.data_collection[spec_1d].world_component_ids
+            wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_first_spec_1d[0]))
 
     app.session.data_collection.add_link(wc_spec_ids)
 
@@ -227,77 +242,59 @@ def mos_spec2d_parser(app, data_obj, data_labels=None, add_to_table=True,
     data_labels : str, optional
         The label applied to the glue data component.
     """
-    # In the case where the data object is a string, attempt to parse it as
-    #  a fits file.
-    # TODO: this current does not handle the case where the file in the path is
-    #  anything but a fits file whose wcs can be extracted.
-    def _parse_as_cube(path):
+    def _parse_as_spectrum1d(path):
+        # Parse as a FITS file and assume the WCS is correct
         with fits.open(path) as hdulist:
             data = hdulist[1].data
             header = hdulist[1].header
-            if header['NAXIS'] == 2:
-                new_data = np.expand_dims(data, axis=1)
-                header['NAXIS'] = 3
-
-            header['NAXIS3'] = 1
-            header['BUNIT'] = 'dN/s'
-            header['CUNIT3'] = 'um'
-            header['CTYPE3'] = 'WAVE'
-
-            # Information not present in the SCI header has to be put there
-            # so spectral_cube won't choke. We cook up a simple linear wcs
-            # with the only intention of making the code run beyond the
-            # spectral_cube processing. There is no guarantee that this will
-            # result in the correct axis label values being displayed.
-            #
-            # This is a stopgap solution that will be replaced when specutils
-            # absorbs the functionality provided by spectral_cube.
-
-            fa = AsdfInFits.open(path)
-            gwcs = fa.tree['meta']['wcs']
-
-            header['CTYPE1'] = 'RA---TAN'
-            header['CTYPE2'] = 'DEC--TAN'
-            header['CUNIT1'] = 'deg'
-            header['CUNIT2'] = 'deg'
-
-            header['CRVAL1'] = gwcs.forward_transform.lon_4.value
-            header['CRVAL2'] = gwcs.forward_transform.lat_4.value
-            header['CRPIX1'] = gwcs.forward_transform.intercept_1.value
-            header['CRPIX2'] = gwcs.forward_transform.intercept_2.value
-            header['CDELT1'] = gwcs.forward_transform.slope_1.value
-            header['CDELT2'] = gwcs.forward_transform.slope_2.value
-            header['PC1_1'] = -1.
-            header['PC1_2'] = 0.
-            header['PC2_1'] = 0.
-            header['PC2_2'] = 1.
-            header['PC3_1'] = 1.
-            header['PC3_2'] = 0.
-
             wcs = WCS(header)
+        return Spectrum1D(data, wcs=wcs)
 
-            meta = {'S_REGION': header['S_REGION'], 'INSTRUME': 'nirspec'}
-
-        return SpectralCube(new_data, wcs=wcs, meta=meta)
-
-    if isinstance(data_labels, str):
-        data_labels = [data_labels]
-
-    # Coerce into list if needed
+    # Coerce into list-like object
     if not isinstance(data_obj, (list, tuple, SpectrumCollection)):
         data_obj = [data_obj]
 
-    data_obj = [_parse_as_cube(x) if _check_is_file(x) else x for x in data_obj]
-
-    if data_labels is None:
-        data_labels = [f"2D Spectrum {i}" for i in range(len(data_obj))]
-    elif len(data_obj) != len(data_labels):
-        data_labels = [f"{data_labels} {i}" for i in range(len(data_obj))]
+    # If we're given a string, repeat it for each object
+    if isinstance(data_labels, str):
+        if len(data_obj) > 1:
+            data_labels = [f"{data_labels} {i}" for i in range(len(data_obj))]
+        else:
+            data_labels = [data_labels]
+    elif data_labels is None:
+        if len(data_obj) > 1:
+            data_labels = [f"2D Spectrum {i}" for i in range(len(data_obj))]
+        else:
+            data_labels = ['2D Spectrum']
 
     with app.data_collection.delay_link_manager_update():
 
-        for i in range(len(data_obj)):
-            app.add_data(data_obj[i], data_labels[i], notify_done=False)
+        for index, data in enumerate(data_obj):
+            # If we got a filepath, first try and parse using the Spectrum1D and
+            # SpectrumList parsers, and then fall back to parsing it as a generic
+            # FITS file.
+            if _check_is_file(data):
+                try:
+                    data = Spectrum1D.read(data)
+                except IORegistryError:
+                    try:
+                        data = Spectrum1D.read(data)
+                    except IORegistryError:
+                        data = _parse_as_spectrum1d(data)
+
+            # Copy (if present) region to top-level meta object
+            if ('header' in data.meta and
+                    'S_REGION' in data.meta['header'] and
+                    'S_REGION' not in data.meta):
+                data.meta['S_REGION'] = data.meta['header']['S_REGION']
+
+            # Set the instrument
+            # TODO: this should not be set to nirspec for all datasets
+            data.meta['INSTRUME'] = 'nirspec'
+
+            # Get the corresponding label for this data product
+            label = data_labels[index]
+
+            app.data_collection[label] = data
 
         if add_to_table:
             _add_to_table(app, data_labels, '2D Spectra')
@@ -542,39 +539,27 @@ def mos_niriss_parser(app, data_dir, obs_label=""):
 
             with fits.open(fname, memmap=False) as temp:
                 sci_hdus = []
+                wav_hdus = {}
                 for i in range(len(temp)):
                     if "EXTNAME" in temp[i].header:
                         if temp[i].header["EXTNAME"] == "SCI":
                             sci_hdus.append(i)
+                            wav_hdus[i] = ('WAVELENGTH', temp[i].header['EXTVER'])
 
-                # Now get a SpectralCube object for each SCI HDU
+                # Now get a Spectrum1D object for each SCI HDU
                 for sci in sci_hdus:
 
                     if temp[sci].header["SPORDER"] == 1:
 
                         data = temp[sci].data
                         meta = temp[sci].header
-                        header = filter_wcs[filter_name]
-
-                        # Information that needs to be added to the header in order to load
-                        # WCS information into SpectralCube.
-                        # TODO: Use gwcs instead to avoid hardcoding information for 3rd wcs axis
-                        new_data = np.expand_dims(data, axis=1)
-                        header['NAXIS'] = 3
-
-                        header['NAXIS3'] = 1
-                        header['BUNIT'] = 'dN/s'
-                        header['CUNIT3'] = 'um'
-                        header['WCSAXES'] = 3
-                        header['CRPIX3'] = 0.0
-                        header['CDELT3'] = 1E-06
-                        header['CUNIT3'] = 'm'
-                        header['CTYPE3'] = 'WAVE'
-                        header['CRVAL3'] = 0.0
-
-                        wcs = WCS(header)
-
-                        spec2d = SpectralCube(new_data, wcs=wcs, meta=meta)
+
+                        # The wavelength is stored in a WAVELENGTH HDU. This is
+                        # a 2D array, but in order to be able to use Spectrum1D
+                        # we use the average wavelength for all image rows
+                        wav = temp[wav_hdus[sci]].data.mean(axis=0) * u.micron
+
+                        spec2d = Spectrum1D(data * u.one, spectral_axis=wav, meta=meta)
 
                         spec2d.meta['INSTRUME'] = 'NIRISS'
 

jdaviz/configs/mosviz/plugins/viewers.py

@@ -3,8 +3,6 @@
 from glue_jupyter.bqplot.profile import BqplotProfileView
 from glue_jupyter.table import TableViewer
 from specutils import Spectrum1D
-from spectral_cube import SpectralCube
-from echo import delay_callback
 import astropy
 from astropy import units as u
 from astropy.utils.introspection import minversion
@@ -78,29 +76,21 @@ def set_plot_axes(self):
 class MosvizProfile2DView(BqplotImageView):
     # Due to limitations in CCDData and 2D data that has spectral and spatial
     #  axes, the default conversion class must handle cubes
-    default_class = SpectralCube
+    default_class = Spectrum1D
 
     tools = ['bqplot:panzoom_x']
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         # Setup viewer option defaults
         self.state.aspect = 'auto'
-        self.state.add_callback('reference_data', self._update_world_axes, priority=100)
 
     def data(self, cls=None):
         return [layer_state.layer.get_object(cls=cls or self.default_class)
                 for layer_state in self.state.layers
                 if hasattr(layer_state, 'layer') and
                 isinstance(layer_state.layer, BaseData)]
 
-    def _update_world_axes(self, data):
-        if data is not None:
-            with delay_callback(self.state, 'x_att_world', 'y_att_world'):
-                if 'Wave' in data.components:
-                    self.state.x_att_world = data.id['Right Ascension']
-                    self.state.y_att_world = data.id['Wave']
-
     def set_plot_axes(self):
         self.figure.axes[0].visible = False
 

jdaviz/configs/mosviz/tests/test_helper.py

@@ -5,7 +5,6 @@
 from astropy.nddata import CCDData
 from astropy.wcs import WCS
 from jdaviz.configs.mosviz.helper import Mosviz
-from spectral_cube import SpectralCube
 from specutils import Spectrum1D, SpectrumCollection
 
 
@@ -27,28 +26,17 @@ def spectrum1d():
 @pytest.fixture
 def spectrum2d():
     header = """
-WCSAXES =                    3 / Number of coordinate axes
-CRPIX1  =               1024.5 / Pixel coordinate of reference point
+WCSAXES =                    2 / Number of coordinate axes
+CRPIX1  =                  0.0 / Pixel coordinate of reference point
 CRPIX2  =               1024.5 / Pixel coordinate of reference point
-CRPIX3  =                  0.0 / Pixel coordinate of reference point
-PC1_1   =                 -1.0 / Coordinate transformation matrix element
-PC3_1   =                  1.0 / Coordinate transformation matrix element
-CDELT1  =  2.8685411111111E-05 / [deg] Coordinate increment at reference point
+CDELT1  =                1E-06 / [m] Coordinate increment at reference point
 CDELT2  =  2.9256727777778E-05 / [deg] Coordinate increment at reference point
-CDELT3  =                1E-06 / [m] Coordinate increment at reference point
-CUNIT1  = 'deg'                / Units of coordinate increment and value
+CUNIT1  = 'm'                  / Units of coordinate increment and value
 CUNIT2  = 'deg'                / Units of coordinate increment and value
-CUNIT3  = 'm'                  / Units of coordinate increment and value
-CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection
-CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection
-CTYPE3  = 'WAVE'               / Vacuum wavelength (linear)
-CRVAL1  =                  5.0 / [deg] Coordinate value at reference point
+CTYPE1  = 'WAVE'               / Vacuum wavelength (linear)
+CTYPE2  = 'OFFSET'             / Spatial offset
+CRVAL1  =                  0.0 / [m] Coordinate value at reference point
 CRVAL2  =                  5.0 / [deg] Coordinate value at reference point
-CRVAL3  =                  0.0 / [m] Coordinate value at reference point
-LONPOLE =                180.0 / [deg] Native longitude of celestial pole
-LATPOLE =                  5.0 / [deg] Native latitude of celestial pole
-MJDREFI =                  0.0 / [d] MJD of fiducial time, integer part
-MJDREFF =                  0.0 / [d] MJD of fiducial time, fractional part
 RADESYS = 'ICRS'               / Equatorial coordinate system
 SPECSYS = 'BARYCENT'           / Reference frame of spectral coordinates
 """
@@ -67,8 +55,8 @@ def spectrum2d():
         new_hdr[key] = value
 
     wcs = WCS(new_hdr)
-    data = np.random.sample((15, 1, 1024))
-    spectral_cube = SpectralCube(data, wcs=wcs)
+    data = np.random.sample((1024, 15)) * u.one
+    spectral_cube = Spectrum1D(data, wcs=wcs)
 
     return spectral_cube
 
@@ -176,6 +164,7 @@ def test_load_list_of_spectrum1d(mosviz_app, spectrum1d):
 
 
 def test_load_spectrum2d(mosviz_app, spectrum2d):
+
     label = "Test 2D Spectrum"
     mosviz_app.load_2d_spectra(spectrum2d, data_labels=label)
 
@@ -187,7 +176,7 @@ def test_load_spectrum2d(mosviz_app, spectrum2d):
 
     data = mosviz_app.app.get_data_from_viewer('spectrum-2d-viewer')
 
-    assert list(data.values())[0].shape == (15, 1, 1024)
+    assert list(data.values())[0].shape == (1024, 15)
     assert list(data.keys())[0] == label
 
 

setup.cfg

@@ -23,7 +23,7 @@ install_requires =
     bqplot>=0.12.29
     bqplot-image-gl>=1.4.1
     glue-core>=1.2.1
-    glue-jupyter>=0.8
+    glue-jupyter>=0.8.1
     echo>=0.5.0
     ipyvue>=1.4.1
     ipyvuetify>=1.7.0
@@ -32,7 +32,7 @@ install_requires =
     voila>=0.2.4
     pyyaml>=5.4.1
     specutils>=1.4.0
-    glue-astronomy>=0.2
+    glue-astronomy>=0.3.2
     click>=7.1.2
     spectral-cube>=0.5
     asteval>=0.9.23