SpectralCube to specutils transition in cubeviz

Concept of loading spectrum1d as cube in CubevizExample

Get slider working correctly

SpectralCube to specutils transition in cubeviz

Manga cube loaded using Spectrum1D, still some things to iron out

Add support for loading Spectrum1D cubes

Fix codestyle checks

Use correct axes without np swapaxes

Set correct RA and Dec to x y axes, respectively

Remove excess commits and get slider working

pllim: Cleaned up logic, remove spectral_cube from parser completely, add tests

adasd

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
@@ -54,10 +54,13 @@ def parse_data(app, file_obj, data_type=None, data_label=None):
     # 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 +70,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 +119,38 @@ 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_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)
 
-    app.add_data(file_obj, data_label)
+        s1d = Spectrum1D(flux=flux, wcs=file_obj.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}")
+        data_label = f"Unknown spectrum object[{attr.upper()}]"
+        app.add_data(s1d, data_label)
 
-    # TODO: SpectralCube does not store mask information
-    # TODO: SpectralCube does not store data quality information
+        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):

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,43 +30,51 @@ 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())
 
     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)))