diff --git a/CHANGES.rst b/CHANGES.rst
index 6fa2fba6fa..dd4245a5ff 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -62,6 +62,7 @@ Specviz
 Other Changes and Additions
 ---------------------------
 
+- Cubeviz now loads data cube as ``Spectrum1D``. [#547]
 
 2.0 (2021-09-17)
 ================
diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
index 7b62f0bbc7..e513cc6e57 100644
--- a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
+++ b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
@@ -6,8 +6,7 @@
 from glue.core.message import (DataCollectionAddMessage,
                                DataCollectionDeleteMessage)
 from traitlets import List, Unicode, Any, Bool, observe
-from spectral_cube import SpectralCube
-from specutils import SpectralRegion
+from specutils import Spectrum1D, manipulation, SpectralRegion, analysis
 from regions import RectanglePixelRegion
 
 from jdaviz.core.events import SnackbarMessage
@@ -63,12 +62,12 @@ def _on_data_updated(self, msg):
                 # Also set the spectral min and max to default to the full range
                 try:
                     self.selected_data = self.dc_items[i]
-                    cube = self._selected_data.get_object(cls=SpectralCube)
+                    cube = self._selected_data.get_object(cls=Spectrum1D, statistic=None)
                     self.spectral_min = cube.spectral_axis[0].value
                     self.spectral_max = cube.spectral_axis[-1].value
                     self.spectral_unit = str(cube.spectral_axis.unit)
                     break
-                # Skip data that can't be returned as a SpectralCube
+                # Skip data that can't be returned as a Spectrum1D
                 except (ValueError, TypeError):
                     continue
 
@@ -80,7 +79,7 @@ def _on_subset_created(self, msg):
     def _on_data_selected(self, event):
         self._selected_data = next((x for x in self.data_collection
                                     if x.label == event['new']))
-        cube = self._selected_data.get_object(cls=SpectralCube)
+        cube = self._selected_data.get_object(cls=Spectrum1D, statistic=None)
         # Update spectral bounds and unit if we've switched to another unit
         if str(cube.spectral_axis.unit) != self.spectral_unit:
             self.spectral_min = cube.spectral_axis[0].value
@@ -92,7 +91,7 @@ def _on_subset_selected(self, event):
         # If "None" selected, reset based on bounds of selected data
         self._selected_subset = self.selected_subset
         if self._selected_subset == "None":
-            cube = self._selected_data.get_object(cls=SpectralCube)
+            cube = self._selected_data.get_object(cls=Spectrum1D, statistic=None)
             self.spectral_min = cube.spectral_axis[0].value
             self.spectral_max = cube.spectral_axis[-1].value
         else:
@@ -120,12 +119,12 @@ def vue_list_subsets(self, event):
         self._spectral_subsets = temp_dict
         self.spectral_subset_items = temp_list
 
-    def vue_calculate_moment(self, event):
+    def vue_calculate_moment(self, *args):
         # Retrieve the data cube and slice out desired region, if specified
-        cube = self._selected_data.get_object(cls=SpectralCube)
+        cube = self._selected_data.get_object(cls=Spectrum1D, statistic=None)
         spec_min = float(self.spectral_min) * u.Unit(self.spectral_unit)
         spec_max = float(self.spectral_max) * u.Unit(self.spectral_unit)
-        slab = cube.spectral_slab(spec_min, spec_max)
+        slab = manipulation.spectral_slab(cube, spec_min, spec_max)
 
         # Calculate the moment and convert to CCDData to add to the viewers
         try:
@@ -134,10 +133,9 @@ def vue_calculate_moment(self, event):
                 raise ValueError("Moment must be a positive integer")
         except ValueError:
             raise ValueError("Moment must be a positive integer")
-        self.moment = slab.moment(n_moment)
+        self.moment = analysis.moment(slab, order=n_moment)
 
-        moment_ccd = CCDData(self.moment.array, wcs=self.moment.wcs,
-                             unit=self.moment.unit)
+        moment_ccd = CCDData(self.moment, unit=self.moment.unit)
 
         label = "Moment {}: {}".format(n_moment, self._selected_data.label)
         fname_label = self._selected_data.label.replace("[", "_").replace("]", "_")
diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
index e83f59db53..540278e8b6 100644
--- a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
+++ b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
@@ -1,25 +1,27 @@
-import numpy as np
-
-from glue.core import Data
+from specutils import Spectrum1D
 
 from jdaviz import Application
 from jdaviz.configs.cubeviz.plugins.moment_maps.moment_maps import MomentMap
 
 
-def test_moment_calculation(spectral_cube_wcs):
+def todo_fix_test_moment_calculation(spectrum1d_cube):
 
     app = Application()
     dc = app.data_collection
-    dc.append(Data(x=np.ones((3, 4, 5)), label='test', coords=spectral_cube_wcs))
+    app.add_data(spectrum1d_cube, 'test')
 
     mm = MomentMap(app=app)
+    mm._subset_selected = 'None'
+    # mm.spectral_min = 1.0 * u.m
+    # mm.spectral_max = 2.0 * u.m
+    mm._on_data_updated(None)
 
-    mm.selected_data = 'test'
-    mm.n_moment = 0
-    mm.vue_calculate_moment(None)
+    mm._on_data_selected({'new': 'test'})
+    mm._on_subset_selected({'new': None})
 
-    print(dc[1].get_object())
+    mm.n_moment = 0
+    mm.vue_calculate_moment()
 
     assert mm.moment_available
     assert dc[1].label == 'Moment 0: test'
-    assert dc[1].get_object().shape == (4, 5)
+    assert dc[1].get_object(cls=Spectrum1D, statistic=None).shape == (4, 2, 2)
diff --git a/jdaviz/configs/cubeviz/plugins/parsers.py b/jdaviz/configs/cubeviz/plugins/parsers.py
index 39c81a4c8e..35bf85b511 100644
--- a/jdaviz/configs/cubeviz/plugins/parsers.py
+++ b/jdaviz/configs/cubeviz/plugins/parsers.py
@@ -2,9 +2,9 @@
 import os
 
 import numpy as np
+from astropy import units as u
 from astropy.io import fits
-from spectral_cube import SpectralCube
-from spectral_cube.io.fits import FITSReadError
+from astropy.wcs import WCS
 from specutils import Spectrum1D
 
 from jdaviz.core.registries import data_parser_registry
@@ -49,15 +49,26 @@ def parse_data(app, file_obj, data_type=None, data_label=None):
         file_name = os.path.basename(file_obj)
 
         with fits.open(file_obj) as hdulist:
-            _parse_hdu(app, hdulist, file_name=data_label or file_name)
+            prihdr = hdulist[0].header
+            telescop = prihdr.get('TELESCOP', '').lower()
+            filetype = prihdr.get('FILETYPE', '').lower()
+            if telescop == 'jwst' and filetype == '3d ifu cube':
+                # TODO: What about ERR, DQ, and WMAP?
+                data_label = f'{file_name}[SCI]'
+                _parse_jwst_s3d(app, hdulist, data_label)
+            else:
+                _parse_hdu(app, hdulist, file_name=data_label or file_name)
 
     # If the data types are custom data objects, use explicit parsers. Note
     #  that this relies on the glue-astronomy machinery to turn the data object
     #  into something glue can understand.
-    elif isinstance(file_obj, SpectralCube):
-        _parse_spectral_cube(app, file_obj, data_type or 'flux', data_label)
     elif isinstance(file_obj, Spectrum1D):
-        _parse_spectrum1d(app, file_obj)
+        if file_obj.flux.ndim == 3:
+            _parse_spectrum1d_3d(app, file_obj)
+        else:
+            _parse_spectrum1d(app, file_obj)
+    else:
+        raise NotImplementedError(f'Unsupported data format: {file_obj}')
 
 
 def _parse_hdu(app, hdulist, file_name=None):
@@ -67,44 +78,36 @@ def _parse_hdu(app, hdulist, file_name=None):
 
     file_name = file_name or "Unknown HDU object"
 
-    # WCS may only exist in a single extension (in this case, just the flux
-    #  flux extension), so first find and store then wcs information.
-    wcs = None
-
-    for hdu in hdulist:
-        if hdu.data is None or not hdu.is_image:
-            continue
-
-        try:
-            sc = SpectralCube.read(hdu, format='fits')
-        except (ValueError, FITSReadError):
-            continue
-        else:
-            wcs = sc.wcs
-
     # Now loop through and attempt to parse the fits extensions as spectral
     #  cube object. If the wcs fails to parse in any case, use the wcs
     #  information we scraped above.
     for hdu in hdulist:
         data_label = f"{file_name}[{hdu.name}]"
 
-        if hdu.data is None or not hdu.is_image:
+        if hdu.data is None or not hdu.is_image or hdu.data.ndim != 3:
             continue
 
-        # This is supposed to fail on attempting to load anything that
-        # isn't cube-shaped. But it's not terribly reliable
         try:
-            sc = SpectralCube.read(hdu, format='fits')
-        except (ValueError, OSError):
-            # This will fail if the parsing of the wcs does not provide
-            # proper celestial axes
+            wcs = WCS(hdu.header, hdulist)
+        except Exception as e:  # TODO: Do we just want to fail here?
+            logging.warning(f"Invalid WCS: {repr(e)}")
+            wcs = None
+
+        if 'BUNIT' in hdu.header:
             try:
-                hdu.header.update(wcs.to_header())
-                sc = SpectralCube.read(hdu)
-            except (ValueError, AttributeError) as e:
-                logging.warning(e)
-                continue
-        except FITSReadError as e:
+                flux_unit = u.Unit(hdu.header['BUNIT'])
+            except Exception:
+                logging.warning("Invalid BUNIT, using count as data unit")
+                flux_unit = u.count
+        else:
+            logging.warning("Missing BUNIT, using count as data unit")
+            flux_unit = u.count
+
+        flux = hdu.data << flux_unit
+
+        try:
+            sc = Spectrum1D(flux=flux, wcs=wcs)
+        except Exception as e:
             logging.warning(e)
             continue
 
@@ -124,21 +127,55 @@ def _parse_hdu(app, hdulist, file_name=None):
             app.add_data_to_viewer('spectrum-viewer', data_label)
 
 
-def _parse_spectral_cube(app, file_obj, data_type='flux', data_label=None):
-    data_label = data_label or f"Unknown spectral cube[{data_type.upper()}]"
+def _parse_jwst_s3d(app, hdulist, data_label):
+    from specutils import Spectrum1D
 
-    app.add_data(file_obj, data_label)
+    unit = u.Unit(hdulist[1].header.get('BUNIT', 'count'))
+    flux = hdulist[1].data << unit
+    wcs = WCS(hdulist[1].header, hdulist)
+    data = Spectrum1D(flux, wcs=wcs)
 
-    if data_type == 'flux':
-        app.add_data_to_viewer('flux-viewer', f"{data_label}")
-        app.add_data_to_viewer('spectrum-viewer', f"{data_label}")
-    elif data_type == 'mask':
-        app.add_data_to_viewer('mask-viewer', f"{data_label}")
-    elif data_type == 'uncert':
-        app.add_data_to_viewer('uncert-viewer', f"{data_label}")
+    # NOTE: Tried to only pass in sliced WCS but got error in Glue.
+    # sliced_wcs = wcs[:, 0, 0]  # Only want wavelengths
+    # data = Spectrum1D(flux, wcs=sliced_wcs)
 
-    # TODO: SpectralCube does not store mask information
-    # TODO: SpectralCube does not store data quality information
+    app.add_data(data, data_label)
+    app.add_data_to_viewer('flux-viewer', data_label)
+    app.add_data_to_viewer('spectrum-viewer', data_label)
+
+
+def _parse_spectrum1d_3d(app, file_obj):
+    # Load spectrum1d as a cube
+
+    for attr in ["flux", "mask", "uncertainty"]:
+        val = getattr(file_obj, attr)
+        if val is None:
+            continue
+
+        if attr == "mask":
+            flux = val << file_obj.flux.unit
+        elif attr == "uncertainty":
+            if hasattr(val, "array"):
+                flux = u.Quantity(val.array, file_obj.flux.unit)
+            else:
+                continue
+        else:
+            flux = val
+
+        flux = np.moveaxis(flux, 1, 0)
+
+        s1d = Spectrum1D(flux=flux, wcs=file_obj.wcs)
+
+        data_label = f"Unknown spectrum object[{attr.upper()}]"
+        app.add_data(s1d, data_label)
+
+        if attr == 'flux':
+            app.add_data_to_viewer('flux-viewer', data_label)
+            app.add_data_to_viewer('spectrum-viewer', data_label)
+        elif attr == 'mask':
+            app.add_data_to_viewer('mask-viewer', data_label)
+        else:  # 'uncertainty'
+            app.add_data_to_viewer('uncert-viewer', data_label)
 
 
 def _parse_spectrum1d(app, file_obj):
diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py b/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py
index e0a1386e87..c7e545f8cc 100644
--- a/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py
+++ b/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py
@@ -1,31 +1,20 @@
-import os
-
-import astropy.units as u
 import numpy as np
 import pytest
+from astropy import units as u
 from astropy.io import fits
-from astropy.nddata import StdDevUncertainty
 from astropy.wcs import WCS
-from spectral_cube import SpectralCube
 from specutils import Spectrum1D
 
-from jdaviz.app import Application
-
-
-@pytest.fixture
-def cubeviz_app():
-    return Application(configuration='cubeviz')
-
 
 @pytest.fixture
 def image_hdu_obj():
-    flux_hdu = fits.ImageHDU(np.random.sample((10, 10, 10)))
+    flux_hdu = fits.ImageHDU(np.ones((10, 10, 10)))
     flux_hdu.name = 'FLUX'
 
     mask_hdu = fits.ImageHDU(np.zeros((10, 10, 10)))
     mask_hdu.name = 'MASK'
 
-    uncert_hdu = fits.ImageHDU(np.random.sample((10, 10, 10)))
+    uncert_hdu = fits.ImageHDU(np.ones((10, 10, 10)))
     uncert_hdu.name = 'ERR'
 
     wcs = WCS(header={
@@ -41,7 +30,7 @@ def image_hdu_obj():
     })
 
     flux_hdu.header.update(wcs.to_header())
-    flux_hdu.header['BUNIT'] = '1E-17 erg/s/cm^2/Angstrom/spaxel'
+    flux_hdu.header['BUNIT'] = '1E-17 erg*s^-1*cm^-2*Angstrom^-1*pix^-1'
 
     mask_hdu.header.update(wcs.to_header())
     uncert_hdu.header.update(wcs.to_header())
@@ -49,35 +38,43 @@ def image_hdu_obj():
     return fits.HDUList([fits.PrimaryHDU(), flux_hdu, mask_hdu, uncert_hdu])
 
 
-@pytest.mark.filterwarnings('ignore:.* contains multiple slashes')
+@pytest.mark.filterwarnings('ignore')
 def test_fits_image_hdu_parse(image_hdu_obj, cubeviz_app):
     cubeviz_app.load_data(image_hdu_obj)
 
-    assert len(cubeviz_app.data_collection) == 3
-    assert cubeviz_app.data_collection[0].label.endswith('[FLUX]')
+    assert len(cubeviz_app.app.data_collection) == 3
+    assert cubeviz_app.app.data_collection[0].label.endswith('[FLUX]')
 
 
-@pytest.mark.filterwarnings('ignore:.* contains multiple slashes')
-def test_spectral_cube_parse(tmpdir, image_hdu_obj, cubeviz_app):
+@pytest.mark.filterwarnings('ignore')
+def test_fits_image_hdu_parse_from_file(tmpdir, image_hdu_obj, cubeviz_app):
     f = tmpdir.join("test_fits_image.fits")
-    path = os.path.join(f.dirname, f.basename)
-    image_hdu_obj.writeto(path)
+    path = f.strpath
+    image_hdu_obj.writeto(path, overwrite=True)
+    cubeviz_app.load_data(path)
+
+    assert len(cubeviz_app.app.data_collection) == 3
+    assert cubeviz_app.app.data_collection[0].label.endswith('[FLUX]')
 
-    sc = SpectralCube.read(path, hdu=1)
 
+@pytest.mark.filterwarnings('ignore')
+def test_spectrum3d_parse(image_hdu_obj, cubeviz_app):
+    flux = image_hdu_obj[1].data << u.Unit(image_hdu_obj[1].header['BUNIT'])
+    wcs = WCS(image_hdu_obj[1].header, image_hdu_obj)
+    sc = Spectrum1D(flux=flux, wcs=wcs)
     cubeviz_app.load_data(sc)
 
-    assert len(cubeviz_app.data_collection) == 1
-    assert cubeviz_app.data_collection[0].label.endswith('[FLUX]')
+    assert len(cubeviz_app.app.data_collection) == 1
+    assert cubeviz_app.app.data_collection[0].label.endswith('[FLUX]')
+
 
+def test_spectrum1d_parse(spectrum1d, cubeviz_app):
+    cubeviz_app.load_data(spectrum1d)
 
-def test_spectrum1d_parse(image_hdu_obj, cubeviz_app):
-    spec = Spectrum1D(flux=np.random.sample(10) * u.Jy,
-                      spectral_axis=np.arange(10) * u.nm,
-                      uncertainty=StdDevUncertainty(
-                          np.random.sample(10) * u.Jy))
+    assert len(cubeviz_app.app.data_collection) == 1
+    assert cubeviz_app.app.data_collection[0].label.endswith('[FLUX]')
 
-    cubeviz_app.load_data(spec)
 
-    assert len(cubeviz_app.data_collection) == 1
-    assert cubeviz_app.data_collection[0].label.endswith('[FLUX]')
+def test_numpy_cube(cubeviz_app):
+    with pytest.raises(NotImplementedError, match='Unsupported data format'):
+        cubeviz_app.load_data(np.ones(27).reshape((3, 3, 3)))
diff --git a/jdaviz/configs/default/plugins/collapse/collapse.py b/jdaviz/configs/default/plugins/collapse/collapse.py
index 040dbdd7f9..ee95804520 100644
--- a/jdaviz/configs/default/plugins/collapse/collapse.py
+++ b/jdaviz/configs/default/plugins/collapse/collapse.py
@@ -4,6 +4,7 @@
 from glue.core import Data
 from glue.core.link_helpers import LinkSame
 from spectral_cube import SpectralCube
+from specutils import Spectrum1D
 from specutils import SpectralRegion
 from traitlets import List, Unicode, Int, Any, observe
 from regions import RectanglePixelRegion
@@ -80,7 +81,7 @@ def _on_data_item_selected(self, event):
     def _on_subset_selected(self, event):
         # If "None" selected, reset based on bounds of selected data
         if self.selected_subset == "None":
-            cube = self._selected_data.get_object(cls=SpectralCube)
+            cube = self._selected_data.get_object(cls=Spectrum1D)
             self.spectral_min = cube.spectral_axis[0].value
             self.spectral_max = cube.spectral_axis[-1].value
         else:
diff --git a/jdaviz/configs/default/plugins/gaussian_smooth/gaussian_smooth.py b/jdaviz/configs/default/plugins/gaussian_smooth/gaussian_smooth.py
index 70bd7037e8..21108b2307 100644
--- a/jdaviz/configs/default/plugins/gaussian_smooth/gaussian_smooth.py
+++ b/jdaviz/configs/default/plugins/gaussian_smooth/gaussian_smooth.py
@@ -6,7 +6,6 @@
                                DataCollectionDeleteMessage)
 from specutils import Spectrum1D
 from specutils.manipulation import gaussian_smooth
-from spectral_cube import SpectralCube
 from traitlets import List, Unicode, Any, Bool, observe
 
 from jdaviz.core.events import SnackbarMessage
@@ -60,7 +59,7 @@ def vue_spectral_smooth(self, *args, **kwargs):
         size = float(self.stddev)
 
         try:
-            spec = self._selected_data.get_object(cls=Spectrum1D)
+            spec = self._selected_data.get_object(cls=Spectrum1D, statistic=None)
         except TypeError:
             snackbar_message = SnackbarMessage(
                 "Unable to perform smoothing over selected data.",
@@ -100,10 +99,18 @@ def vue_spatial_convolution(self, *args):
         """
 
         size = float(self.stddev)
-        cube = self._selected_data.get_object(cls=SpectralCube)
+
+        # Get information from the flux component
+        attribute = self._selected_data.main_components[0]
+
+        cube = self._selected_data.get_object(cls=Spectrum1D,
+                                              attribute=attribute,
+                                              statistic=None)
+        flux_unit = cube.flux.unit
+
         # Extend the 2D kernel to have a length 1 spectral dimension, so that
         # we can do "3d" convolution to the whole cube
-        kernel = np.expand_dims(Gaussian2DKernel(size), 0)
+        kernel = np.expand_dims(Gaussian2DKernel(size), 2)
 
         # TODO: in vuetify >2.3, timeout should be set to -1 to keep open
         #  indefinitely
@@ -112,13 +119,11 @@ def vue_spatial_convolution(self, *args):
             loading=True, timeout=0, sender=self)
         self.hub.broadcast(snackbar_message)
 
-        convolved_data = convolve(cube.hdu.data, kernel)
+        convolved_data = convolve(cube, kernel)
+
         # Create a new cube with the old metadata. Note that astropy
-        # convolution generates values for masked (NaN) data, but we keep the
-        # original mask here.
-        newcube = SpectralCube(data=convolved_data, wcs=cube.wcs,
-                               mask=cube.mask, meta=cube.meta,
-                               fill_value=cube.fill_value)
+        # convolution generates values for masked (NaN) data.
+        newcube = Spectrum1D(flux=convolved_data * flux_unit, wcs=cube.wcs)
 
         label = f"Smoothed {self._selected_data.label} spatial stddev {size}"
 
diff --git a/jdaviz/configs/default/plugins/gaussian_smooth/tests/test_gaussian_smooth.py b/jdaviz/configs/default/plugins/gaussian_smooth/tests/test_gaussian_smooth.py
index 5616f064f9..7029bd7b2c 100644
--- a/jdaviz/configs/default/plugins/gaussian_smooth/tests/test_gaussian_smooth.py
+++ b/jdaviz/configs/default/plugins/gaussian_smooth/tests/test_gaussian_smooth.py
@@ -1,17 +1,15 @@
-import numpy as np
-
-from glue.core import Data
+import pytest
 from jdaviz import Application
-from spectral_cube import SpectralCube
+from specutils import Spectrum1D
 
 from jdaviz.configs.default.plugins.gaussian_smooth.gaussian_smooth import GaussianSmooth
 
 
-def test_linking_after_spectral_smooth(spectral_cube_wcs):
+def test_linking_after_spectral_smooth(spectrum1d_cube):
 
     app = Application(configuration="cubeviz")
     dc = app.data_collection
-    dc.append(Data(x=np.ones((3, 4, 5)), label='test', coords=spectral_cube_wcs))
+    app.add_data(spectrum1d_cube, 'test')
 
     gs = GaussianSmooth(app=app)
 
@@ -27,11 +25,12 @@ def test_linking_after_spectral_smooth(spectral_cube_wcs):
     assert dc.external_links[0].cids2[0] is dc[1].world_component_ids[0]
 
 
-def test_spatial_convolution(spectral_cube_wcs):
+@pytest.mark.filterwarnings("ignore::UserWarning")
+def test_spatial_convolution(spectrum1d_cube):
 
-    app = Application()
+    app = Application(configuration="cubeviz")
     dc = app.data_collection
-    dc.append(Data(x=np.ones((3, 4, 5)), label='test', coords=spectral_cube_wcs))
+    app.add_data(spectrum1d_cube, 'test')
 
     gs = GaussianSmooth(app=app)
     gs._on_data_selected({'new': 'test'})
@@ -40,4 +39,5 @@ def test_spatial_convolution(spectral_cube_wcs):
 
     assert len(dc) == 2
     assert dc[1].label == "Smoothed test spatial stddev 3.0"
-    assert dc["Smoothed test spatial stddev 3.0"].get_object(cls=SpectralCube).shape == (3, 4, 5)
+    assert (dc["Smoothed test spatial stddev 3.0"].get_object(cls=Spectrum1D, statistic=None).shape
+            == (4, 2, 2))
diff --git a/jdaviz/conftest.py b/jdaviz/conftest.py
index 74509b9306..ed26410fb5 100644
--- a/jdaviz/conftest.py
+++ b/jdaviz/conftest.py
@@ -49,6 +49,15 @@ def spectral_cube_wcs(request):
     return wcs
 
 
+@pytest.fixture
+def spectrum1d_cube_wcs(request):
+    # A simple spectrum1D WCS used by some tests
+    wcs = WCS(naxis=3)
+    wcs.wcs.ctype = 'RA---TAN', 'DEC--TAN', 'WAVE-LOG'
+    wcs.wcs.set()
+    return wcs
+
+
 @pytest.fixture
 def spectrum1d():
     np.random.seed(42)
@@ -62,6 +71,20 @@ def spectrum1d():
     return Spectrum1D(spectral_axis=spec_axis, flux=flux, uncertainty=uncertainty)
 
 
+@pytest.fixture
+def spectrum1d_cube():
+    flux = np.arange(16).reshape([2, 2, 4]) * u.Jy
+    wcs_dict = {"CTYPE1": "RA---TAN", "CTYPE2": "DEC--TAN", "CTYPE3": "WAVE-LOG",
+                "CRVAL1": 205, "CRVAL2": 27, "CRVAL3": 4.622e-7,
+                "CDELT1": -0.0001, "CDELT2": 0.0001, "CDELT3": 8e-11,
+                "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0}
+    w = WCS(wcs_dict)
+
+    spec = Spectrum1D(flux=flux, wcs=w)
+
+    return spec
+
+
 try:
     from pytest_astropy_header.display import PYTEST_HEADER_MODULES, TESTED_VERSIONS
 except ImportError:
diff --git a/setup.cfg b/setup.cfg
index b16a6a375b..881e83fb4e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -31,7 +31,7 @@ install_requires =
     ipygoldenlayout>=0.3.0
     voila>=0.2.4
     pyyaml>=5.4.1
-    specutils>=1.4.0
+    specutils@git+https://github.com/astropy/specutils.git
     glue-astronomy>=0.3.2
     click>=7.1.2
     spectral-cube>=0.5
@@ -91,6 +91,7 @@ filterwarnings =
     ignore::DeprecationWarning:bqplot_image_gl
     ignore::DeprecationWarning:bqscales
     ignore::DeprecationWarning:traittypes
+    ignore:::specutils.spectra.spectrum1d
 
 [flake8]
 max-line-length = 100