toggle flux <-> surface brightness units (#2781)

* add toggle feature between surface brightness and flux, currently incomplete

* custom eqv added, hide toggle with traitlet, move test

* add sb units to equiv_units(), cursor to recognize default data toggle, sb/flux dropdown updating with correct units

* add sb units to valid units, change from toggle to dropdown for sb/flux

* changing from switch to dropdown, updating logic for conversion

* removing toggle from spectral extraction

* resolve UnitConversionError tracebacks in notebook

* change flux_unit to flux_or_sb_unit in tests, seperate translation test from sb conversion test

* resolve CI test failures

* adding flux_or_sb doc entry, and TEMP exposing flux_or_sb_unit to verify CI testing

* fix styling

* add test coverage

* add change log

CHANGES.rst

@@ -4,6 +4,9 @@
 New Features
 ------------
 
+- Adding flux/surface brightness translation and surface brightness
+  unit conversion in Cubeviz and Specviz. [#2781]
+
 Cubeviz
 ^^^^^^^
 

jdaviz/app.py

@@ -66,20 +66,29 @@
 
 @unit_converter('custom-jdaviz')
 class UnitConverterWithSpectral:
-
     def equivalent_units(self, data, cid, units):
         if cid.label == "flux":
             eqv = u.spectral_density(1 * u.m)  # Value does not matter here.
             list_of_units = set(list(map(str, u.Unit(units).find_equivalent_units(
-                include_prefix_units=True, equivalencies=eqv))) + [
-                    'Jy', 'mJy', 'uJy',
+                include_prefix_units=True, equivalencies=eqv)))
+                + [
+                    'Jy', 'mJy', 'uJy', 'MJy',
                     'W / (m2 Hz)', 'W / (Hz m2)',  # Order is different in astropy v5.3
                     'eV / (s m2 Hz)', 'eV / (Hz s m2)',
                     'erg / (s cm2)',
-                    'erg / (s cm2 Angstrom)', 'erg / (Angstrom s cm2)',
+                    'erg / (s cm2 Angstrom)', 'erg / (s cm2 Angstrom)',
                     'erg / (s cm2 Hz)', 'erg / (Hz s cm2)',
-                    'ph / (s cm2 Angstrom)', 'ph / (Angstrom s cm2)',
-                    'ph / (s cm2 Hz)', 'ph / (Hz s cm2)'
+                    'ph / (s cm2 Angstrom)', 'ph / (s cm2 Angstrom)',
+                    'ph / (Hz s cm2)', 'ph / (Hz s cm2)', 'bol', 'AB', 'ST'
+                ]
+                + [
+                    'Jy / sr', 'mJy / sr', 'uJy / sr', 'MJy / sr',
+                    'W / (Hz sr m2)',
+                    'eV / (s m2 Hz sr)',
+                    'erg / (s cm2 sr)',
+                    'erg / (s cm2 Angstrom sr)', 'erg / (s cm2 Hz sr)',
+                    'ph / (s cm2 Angstrom sr)', 'ph / (s cm2 Hz sr)',
+                    'bol / sr', 'AB / sr', 'ST / sr'
                 ])
         else:  # spectral axis
             # prefer Hz over Bq and um over micron
@@ -100,12 +109,48 @@ def to_unit(self, data, cid, values, original_units, target_units):
             except RuntimeError:
                 eqv = []
             else:
-                if len(values) == 2:
+                # Ensure a spectrum passed through Spectral Extraction plugin
+                if '_pixel_scale_factor' in spec.meta:
+                    # if spectrum data collection item is in Surface Brightness units
+                    if u.sr in spec.unit.bases:
+                        # Data item in data collection does not update from conversion/translation.
+                        # App wide orginal data units are used for conversion, orginal_units and
+                        # target_units dicate the conversion to take place.
+                        if (u.sr in u.Unit(original_units).bases) and \
+                           (u.sr not in u.Unit(target_units).bases):
+                            # Surface Brightness -> Flux
+                            eqv = [(u.MJy / u.sr,
+                                    u.MJy,
+                                    lambda x: (x * spec.meta['_pixel_scale_factor']),
+                                    lambda x: x)]
+                        else:
+                            # Flux -> Surface Brightness
+                            eqv = u.spectral_density(spec.spectral_axis)
+
+                    # if spectrum data collection item is in Flux units
+                    elif u.sr not in spec.unit.bases:
+                        # Data item in data collection does not update from conversion/translation.
+                        # App wide orginal data units are used for conversion, orginal_units and
+                        # target_units dicate the conversion to take place.
+                        if (u.sr not in u.Unit(original_units).bases) and \
+                           (u.sr in u.Unit(target_units).bases):
+                            # Flux -> Surface Brightness
+                            eqv = [(u.MJy,
+                                    u.MJy / u.sr,
+                                    lambda x: (x / spec.meta['_pixel_scale_factor']),
+                                    lambda x: x)]
+                        else:
+                            # Surface Brightness -> Flux
+                            eqv = u.spectral_density(spec.spectral_axis)
+
+                elif len(values) == 2:
                     # Need this for setting the y-limits
                     spec_limits = [spec.spectral_axis[0].value, spec.spectral_axis[-1].value]
                     eqv = u.spectral_density(spec_limits * spec.spectral_axis.unit)
+
                 else:
                     eqv = u.spectral_density(spec.spectral_axis)
+
         else:  # spectral axis
             eqv = u.spectral()
 

jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py

@@ -3,7 +3,6 @@
 
 import numpy as np
 import astropy
-from astropy import units as u
 from astropy.utils.decorators import deprecated
 from astropy.nddata import (
     NDDataArray, StdDevUncertainty
@@ -562,13 +561,3 @@ def _live_update(self, event={}):
         for mark in self.marks.values():
             mark.update_xy(sp.spectral_axis.value, sp.flux.value)
             mark.visible = True
-
-    def translate_units(self, collapsed_spec):
-        # remove sr
-        if u.sr in collapsed_spec._unit.bases:
-            collapsed_spec._data *= collapsed_spec.meta['_pixel_scale_factor']
-            collapsed_spec._unit *= u.sr
-        # add sr
-        elif u.sr not in collapsed_spec._unit.bases:
-            collapsed_spec._data /= collapsed_spec.meta['_pixel_scale_factor']
-            collapsed_spec._unit /= u.sr

jdaviz/configs/cubeviz/plugins/spectral_extraction/tests/test_spectral_extraction.py

@@ -12,7 +12,6 @@
 from regions import (CirclePixelRegion, CircleAnnulusPixelRegion, EllipsePixelRegion,
                      RectanglePixelRegion, PixCoord)
 from specutils import Spectrum1D
-from astropy.wcs import WCS
 
 
 def test_version_after_nddata_update(cubeviz_helper, spectrum1d_cube_with_uncerts):
@@ -375,53 +374,6 @@ def test_cube_extraction_with_nan(cubeviz_helper, image_cube_hdu_obj):
     assert_allclose(sp_subset.flux.value, 12)  # (4 x 4) - 4
 
 
-def test_unit_translation(cubeviz_helper):
-    # custom cube so we have PIXAR_SR in metadata, and flux units = Jy/pix
-    wcs_dict = {"CTYPE1": "WAVE-LOG", "CTYPE2": "DEC--TAN", "CTYPE3": "RA---TAN",
-                "CRVAL1": 4.622e-7, "CRVAL2": 27, "CRVAL3": 205,
-                "CDELT1": 8e-11, "CDELT2": 0.0001, "CDELT3": -0.0001,
-                "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0, "PIXAR_SR": 8e-11}
-    w = WCS(wcs_dict)
-    flux = np.zeros((30, 20, 3001), dtype=np.float32)
-    flux[5:15, 1:11, :] = 1
-    cube = Spectrum1D(flux=flux * u.MJy, wcs=w, meta=wcs_dict)
-    cubeviz_helper.load_data(cube, data_label="test")
-
-    center = PixCoord(5, 10)
-    cubeviz_helper.load_regions(CirclePixelRegion(center, radius=2.5))
-
-    extract_plg = cubeviz_helper.plugins['Spectral Extraction']
-
-    extract_plg.aperture = extract_plg.aperture.choices[-1]
-    extract_plg.aperture_method.selected = 'Exact'
-    extract_plg.wavelength_dependent = True
-    extract_plg.function = 'Sum'
-    # set so pixel scale factor != 1
-    extract_plg.reference_spectral_value = 0.000001
-
-    # collapse to spectrum, now we can get pixel scale factor
-    collapsed_spec = extract_plg.collapse_to_spectrum()
-
-    assert collapsed_spec.meta['_pixel_scale_factor'] != 1
-
-    # store to test second time after calling translate_units
-    mjy_sr_data1 = collapsed_spec._data[0]
-
-    extract_plg._obj.translate_units(collapsed_spec)
-
-    assert collapsed_spec._unit == u.MJy / u.sr
-    # some value in MJy/sr that we know the outcome after translation
-    assert np.allclose(collapsed_spec._data[0], 8.7516529e10)
-
-    extract_plg._obj.translate_units(collapsed_spec)
-
-    # translating again returns the original units
-    assert collapsed_spec._unit == u.MJy
-    # returns to the original values
-    # which is a value in Jy/pix that we know the outcome after translation
-    assert np.allclose(collapsed_spec._data[0], mjy_sr_data1)
-
-
 def test_autoupdate_results(cubeviz_helper, spectrum1d_cube_largest):
     cubeviz_helper.load_data(spectrum1d_cube_largest)
     fv = cubeviz_helper.viewers['flux-viewer']._obj

jdaviz/configs/imviz/plugins/coords_info/coords_info.py

@@ -573,8 +573,18 @@ def _copy_axes_to_spectral():
 
                 # Calculations have to happen in the frame of viewer display units.
                 disp_wave = sp.spectral_axis.to_value(viewer.state.x_display_unit, u.spectral())
-                disp_flux = sp.flux.to_value(viewer.state.y_display_unit,
-                                             u.spectral_density(sp.spectral_axis))
+
+                # temporarily here, may be removed after upstream units handling
+                # or will be generalized for any sb <-> flux
+                if '_pixel_scale_factor' in sp.meta:
+                    eqv = [(u.MJy / u.sr,
+                            u.MJy,
+                            lambda x: (x * sp.meta['_pixel_scale_factor']),
+                            lambda x: x)]
+                    disp_flux = sp.flux.to_value(viewer.state.y_display_unit, eqv)
+                else:
+                    disp_flux = sp.flux.to_value(viewer.state.y_display_unit,
+                                                 u.spectral_density(sp.spectral_axis))
 
                 # Out of range in spectral axis.
                 if (self.dataset.selected != lyr.layer.label and

jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py

@@ -5,6 +5,8 @@
 from astropy.nddata import InverseVariance
 from specutils import Spectrum1D
 from astropy.utils.introspection import minversion
+from astropy.wcs import WCS
+from regions import PixCoord, CirclePixelRegion
 
 ASTROPY_LT_5_3 = not minversion(astropy, "5.3")
 
@@ -120,3 +122,113 @@ def test_non_stddev_uncertainty(specviz_helper):
         np.abs(viewer.figure.marks[-1].y - viewer.figure.marks[-1].y.mean(0)),
         stddev
     )
+
+
+def test_unit_translation(cubeviz_helper):
+    # custom cube so PIXAR_SR is in metadata, and Flux units, and in MJy
+    wcs_dict = {"CTYPE1": "WAVE-LOG", "CTYPE2": "DEC--TAN", "CTYPE3": "RA---TAN",
+                "CRVAL1": 4.622e-7, "CRVAL2": 27, "CRVAL3": 205,
+                "CDELT1": 8e-11, "CDELT2": 0.0001, "CDELT3": -0.0001,
+                "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0, "PIXAR_SR": 8e-11}
+    w = WCS(wcs_dict)
+    flux = np.zeros((30, 20, 3001), dtype=np.float32)
+    flux[5:15, 1:11, :] = 1
+    cube = Spectrum1D(flux=flux * u.MJy, wcs=w, meta=wcs_dict)
+    cubeviz_helper.load_data(cube, data_label="test")
+
+    center = PixCoord(5, 10)
+    cubeviz_helper.load_regions(CirclePixelRegion(center, radius=2.5))
+
+    uc_plg = cubeviz_helper.plugins['Unit Conversion']
+    # we can get rid of this after all spectra pass through
+    # spectral extraction plugin
+    extract_plg = cubeviz_helper.plugins['Spectral Extraction']
+
+    extract_plg.aperture = extract_plg.aperture.choices[-1]
+    extract_plg.aperture_method.selected = 'Exact'
+    extract_plg.wavelength_dependent = True
+    extract_plg.function = 'Sum'
+    # set so pixel scale factor != 1
+    extract_plg.reference_spectral_value = 0.000001
+
+    # all spectra will pass through spectral extraction,
+    # this will store a scale factor for use in translations.
+    collapsed_spec = extract_plg.collapse_to_spectrum()
+
+    # test that the scale factor was set
+    assert collapsed_spec.meta['_pixel_scale_factor'] != 1
+
+    # When the dropdown is displayed, this ensures the loaded
+    # data collection item units will be used for translations.
+    uc_plg._obj.show_translator = True
+    assert uc_plg._obj.flux_or_sb_selected == 'Flux'
+
+    # to have access to display units
+    viewer_1d = cubeviz_helper.app.get_viewer(
+        cubeviz_helper._default_spectrum_viewer_reference_name)
+
+    #  for testing _set_flux_or_sb()
+    uc_plg._obj.show_translator = False
+
+    # change global y-units from Flux -> Surface Brightness
+    uc_plg._obj.flux_or_sb_selected = 'Surface Brightness'
+
+    uc_plg._obj.show_translator = True
+    assert uc_plg._obj.flux_or_sb_selected == 'Surface Brightness'
+    y_display_unit = u.Unit(viewer_1d.state.y_display_unit)
+
+    # check if units translated
+    assert y_display_unit == u.MJy / u.sr
+
+
+def test_sb_unit_conversion(cubeviz_helper):
+    # custom cube to have Surface Brightness units
+    wcs_dict = {"CTYPE1": "WAVE-LOG", "CTYPE2": "DEC--TAN", "CTYPE3": "RA---TAN",
+                "CRVAL1": 4.622e-7, "CRVAL2": 27, "CRVAL3": 205,
+                "CDELT1": 8e-11, "CDELT2": 0.0001, "CDELT3": -0.0001,
+                "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0, "PIXAR_SR": 8e-11}
+    w = WCS(wcs_dict)
+    flux = np.zeros((30, 20, 3001), dtype=np.float32)
+    flux[5:15, 1:11, :] = 1
+    cube = Spectrum1D(flux=flux * (u.MJy / u.sr), wcs=w, meta=wcs_dict)
+    cubeviz_helper.load_data(cube, data_label="test")
+
+    uc_plg = cubeviz_helper.plugins['Unit Conversion']
+    uc_plg.open_in_tray()
+
+    # to have access to display units
+    viewer_1d = cubeviz_helper.app.get_viewer(
+        cubeviz_helper._default_spectrum_viewer_reference_name)
+
+    # Surface Brightness conversion
+    uc_plg.flux_or_sb_unit = 'Jy / sr'
+    y_display_unit = u.Unit(viewer_1d.state.y_display_unit)
+    assert y_display_unit == u.Jy / u.sr
+
+    # Try a second conversion
+    uc_plg.flux_or_sb_unit = 'W / Hz sr m2'
+    y_display_unit = u.Unit(viewer_1d.state.y_display_unit)
+    assert y_display_unit == u.Unit("W / (Hz sr m2)")
+
+    # really a translation test, test_unit_translation loads a Flux
+    # cube, this test load a Surface Brightness Cube, this ensures
+    # two-way translation
+    uc_plg.flux_or_sb_unit = 'MJy / sr'
+    y_display_unit = u.Unit(viewer_1d.state.y_display_unit)
+
+    # we can get rid of this after all spectra pass through
+    # spectral extraction plugin
+    extract_plg = cubeviz_helper.plugins['Spectral Extraction']
+    extract_plg.aperture = extract_plg.aperture.choices[-1]
+    extract_plg.aperture_method.selected = 'Exact'
+    extract_plg.wavelength_dependent = True
+    extract_plg.function = 'Sum'
+    extract_plg.reference_spectral_value = 0.000001
+    extract_plg.collapse_to_spectrum()
+
+    uc_plg._obj.show_translator = True
+    uc_plg._obj.flux_or_sb_selected = 'Flux'
+    uc_plg.flux_or_sb_unit = 'MJy'
+    y_display_unit = u.Unit(viewer_1d.state.y_display_unit)
+
+    assert y_display_unit == u.MJy