Correct sum units for per steradian (spacetelescope#2873)

* Correct sum units for per steradian

* Rebase with handle_display_unit fix and update tests

* Fix moment map test

* Fix remote data test

* Change _flux and _spectral conversions to be app level

* Fix style

* Move flux and spectral axis conversions to utils

* Remove unused import

* Remove other unused import

* Add description of flux conversion method

* Update changes file

* Update changes file again

CHANGES.rst

@@ -96,6 +96,9 @@ Cubeviz
 
 - Handle display units when doing flux / surface brightness conversions. [#2910]
 
+- Flux units are now correct for collapsed spectra when using the sum function
+  when units are in per steradian. [#2873]
+
 Imviz
 ^^^^^
 

jdaviz/app.py

@@ -9,7 +9,6 @@
 from astropy.nddata import NDData, NDDataArray
 from astropy.io import fits
 from astropy.time import Time
-from astropy.units import Quantity
 from echo import CallbackProperty, DictCallbackProperty, ListCallbackProperty
 from ipygoldenlayout import GoldenLayout
 from ipysplitpanes import SplitPanes
@@ -51,7 +50,8 @@
                                     data_parser_registry)
 from jdaviz.core.tools import ICON_DIR
 from jdaviz.utils import (SnackbarQueue, alpha_index, data_has_valid_wcs, layer_is_table_data,
-                          MultiMaskSubsetState, _wcs_only_label)
+                          MultiMaskSubsetState, _wcs_only_label, flux_conversion,
+                          spectral_axis_conversion)
 
 __all__ = ['Application', 'ALL_JDAVIZ_CONFIGS', 'UnitConverterWithSpectral']
 
@@ -111,74 +111,9 @@ def to_unit(self, data, cid, values, original_units, target_units):
             except RuntimeError:
                 data = data.get_object(cls=NDDataArray)
                 spec = Spectrum1D(flux=data.data * u.Unit(original_units))
-            return self._flux_conversion(spec, values, original_units, target_units)
+            return flux_conversion(spec, values, original_units, target_units)
         else:  # spectral axis
-            return self._spectral_axis_conversion(values, original_units, target_units)
-
-    def _flux_conversion(self, spec, values, original_units, target_units):
-        # If there are only two values, this is likely the limits being converted, so then
-        # in case we need to use the spectral density equivalency, we need to provide only
-        # to spectral axis values. If there is only one value
-        if not np.isscalar(values) and len(values) == 2:
-            spectral_values = spec.spectral_axis[0]
-        else:
-            spectral_values = spec.spectral_axis
-
-        # Ensure a spectrum passed through Spectral Extraction plugin
-        if '_pixel_scale_factor' in spec.meta and len(values) != 2:
-            # Data item in data collection does not update from conversion/translation.
-            # App wide original data units are used for conversion, original and
-            # target_units dictate 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 * np.array(spec.meta.get('_pixel_scale_factor', 1))),
-                        lambda x: x)]
-            elif (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 / np.array(spec.meta.get('_pixel_scale_factor', 1))),
-                        lambda x: x)]
-            else:
-                eqv = u.spectral_density(spectral_values)
-
-        elif len(values) == 2:
-            # Need this for setting the y-limits
-            eqv = u.spectral_density(spectral_values)
-
-            if '_pixel_scale_factor' in spec.meta:
-                # get min and max scale factors, to use with min and max of spec for
-                # y-limits. Make sure they are Quantities (can be numpy.float64).
-                pixel_scale_min = (Quantity(min(spec.meta.get('_pixel_scale_factor', 1)))).value
-                pixel_scale_max = (Quantity(max(spec.meta.get('_pixel_scale_factor', 1)))).value
-                min_max = [pixel_scale_min, pixel_scale_max]
-
-                if (u.sr in u.Unit(original_units).bases) and \
-                        (u.sr not in u.Unit(target_units).bases):
-                    eqv += [(u.MJy,
-                             u.MJy / u.sr,
-                             lambda x: x * np.array(min_max),
-                             lambda x: x)]
-                elif (u.sr not in u.Unit(original_units).bases) and \
-                        (u.sr in u.Unit(target_units).bases):
-                    eqv += [(u.MJy / u.sr,
-                             u.MJy,
-                             lambda x: x / np.array(min_max),
-                             lambda x: x)]
-
-        else:
-            eqv = u.spectral_density(spectral_values)
-
-        return (values * u.Unit(original_units)).to_value(u.Unit(target_units), equivalencies=eqv)
-
-    def _spectral_axis_conversion(self, values, original_units, target_units):
-        eqv = u.spectral() + u.pixel_scale(1*u.pix)
-        return (values * u.Unit(original_units)).to_value(u.Unit(target_units), equivalencies=eqv)
+            return spectral_axis_conversion(values, original_units, target_units)
 
 
 # Set default opacity for data layers to 1 instead of 0.8 in

jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py

@@ -22,7 +22,7 @@
                                         with_spinner)
 from jdaviz.core.validunits import check_if_unit_is_per_solid_angle
 from jdaviz.core.user_api import PluginUserApi
-from jdaviz.app import UnitConverterWithSpectral as uc
+from jdaviz.utils import flux_conversion
 
 
 __all__ = ['MomentMap']
@@ -366,9 +366,9 @@ def calculate_moment(self, add_data=True):
             pix_scale = self.dataset.selected_dc_item.meta.get('PIXAR_SR', 1.0)
             pix_scale_factor = (flux_values * pix_scale)
             temp_spec.meta['_pixel_scale_factor'] = pix_scale_factor
-            converted_spec = uc._flux_conversion(self, temp_spec, self.moment.value,
-                                                 self.moment.unit,
-                                                 moment_new_unit) * moment_new_unit
+            converted_spec = flux_conversion(temp_spec, self.moment.value,
+                                             self.moment.unit,
+                                             moment_new_unit) * moment_new_unit
             self.moment = converted_spec
 
         # Reattach the WCS so we can load the result

jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py

@@ -277,13 +277,19 @@ def test_momentmap_nirspec_prism(cubeviz_helper, tmp_path):
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
         cubeviz_helper.load_data(uri, cache=True, local_path=local_path)
-    plugin = cubeviz_helper.plugins['Moment Maps']
-    plugin.calculate_moment()
-    assert isinstance(plugin._obj.moment.wcs, WCS)
+    uc = cubeviz_helper.plugins["Unit Conversion"]
+    uc.open_in_tray()  # plugin has to be open for unit change to take hold
+    uc._obj.show_translator = True
+    uc.flux_or_sb.selected = 'Surface Brightness'
+    mm = cubeviz_helper.plugins['Moment Maps']._obj
+    mm.open_in_tray()  # plugin has to be open for unit change to take hold
+    mm._set_data_units()
+    mm.calculate_moment()
+    assert isinstance(mm.moment.wcs, WCS)
 
     # Because cube axes order is re-arranged by specutils on load, this gets confusing.
     # There is a swapaxes within Moment Map WCS calculation.
-    sky_moment = plugin._obj.moment.wcs.pixel_to_world(50, 30)
+    sky_moment = mm.moment.wcs.pixel_to_world(50, 30)
     sky_cube = cubeviz_helper.app.data_collection[0].meta["_orig_spec"].wcs.celestial.pixel_to_world(30, 50)  # noqa: E501
     assert_allclose((sky_moment.ra.deg, sky_moment.dec.deg),
                     (sky_cube.ra.deg, sky_cube.dec.deg))
@@ -307,8 +313,13 @@ def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_
         warnings.filterwarnings("ignore", message="No observer defined on WCS.*")
         cubeviz_helper.load_data(sb_cube, data_label='test')
 
+    uc = cubeviz_helper.plugins["Unit Conversion"]
+    uc.open_in_tray()  # plugin has to be open for unit change to take hold
+    uc._obj.show_translator = True
+    uc.flux_or_sb.selected = 'Surface Brightness'
     mm = cubeviz_helper.plugins['Moment Maps']._obj
     mm.open_in_tray()  # plugin has to be open for unit change to take hold
+    mm._set_data_units()
 
     # check that label is initialized with 'Surface Brightness' since the cube
     # loaded is in MJy / sr. for the 0th moment, the only item will be the 0th
@@ -322,8 +333,7 @@ def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_
     assert mm.moment.unit == f'M{moment_unit}'
 
     # now change surface brightness units in the unit conversion plugin
-    uc = cubeviz_helper.plugins["Unit Conversion"]
-    uc.open_in_tray()  # plugin has to be open for unit change to take hold
+
     uc.flux_or_sb_unit = 'Jy / sr'
 
     # and make sure this change is propogated
@@ -335,7 +345,6 @@ def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_
     mm.calculate_moment()
     assert mm.moment.unit == moment_unit
 
-    uc._obj.show_translator = True
     uc.flux_or_sb.selected = 'Flux'
     mm._set_data_units()
 

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

@@ -3,6 +3,7 @@
 
 import numpy as np
 import astropy
+import astropy.units as u
 from astropy.nddata import (
     NDDataArray, StdDevUncertainty
 )
@@ -426,11 +427,20 @@ def _extract_from_aperture(self, spectral_cube, uncert_cube, aperture,
                                            uncertainty=collapsed_as_mean.uncertainty,
                                            wcs=collapsed_as_mean.wcs,
                                            meta=collapsed_as_mean.meta)
+        elif selected_func == 'sum':
+            collapsed_nddata = getattr(nddata_reshaped, selected_func)(
+                axis=spatial_axes, **kwargs
+            )  # returns an NDDataArray
+            # Remove per steradian denominator
+            if astropy.units.sr in collapsed_nddata.unit.bases:
+                aperture_area = (self.aperture_area_along_spectral
+                                 * self.spectral_cube.meta.get('PIXAR_SR', 1.0) * u.sr)
+                collapsed_nddata = collapsed_nddata.multiply(aperture_area,
+                                                             propagate_uncertainties=True)
         else:
             collapsed_nddata = getattr(nddata_reshaped, selected_func)(
                 axis=spatial_axes, **kwargs
             )  # returns an NDDataArray
-
         # Convert to Spectrum1D, with the spectral axis in correct units:
         if hasattr(spectral_cube.coords, 'spectral_wcs'):
             target_wave_unit = spectral_cube.coords.spectral_wcs.world_axis_units[0]

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

@@ -504,3 +504,17 @@ def test_extraction_composite_subset(cubeviz_helper, spectrum1d_cube):
         (spectrum_1 + spectrum_2).flux.value,
         (spectrum_3).flux.value
     )
+
+
+def test_spectral_extraction_with_correct_sum_units(cubeviz_helper,
+                                                    spectrum1d_cube_fluxunit_jy_per_steradian):
+    cubeviz_helper.load_data(spectrum1d_cube_fluxunit_jy_per_steradian)
+    spec_extr_plugin = cubeviz_helper.plugins['Spectral Extraction']._obj
+    collapsed = spec_extr_plugin.extract()
+    np.testing.assert_allclose(
+        collapsed.flux.value,
+        [3800., 11800., 19800., 27800., 35800., 43800., 51800., 59800., 67800., 75800.]
+
+    )
+    assert collapsed.flux.unit == u.Jy
+    assert collapsed.uncertainty.unit == u.Jy

jdaviz/configs/specviz/plugins/line_analysis/tests/test_lineflux.py

@@ -12,9 +12,9 @@
 
 # Maps input spectrum flux unit to expected line analysis flux unit
 expected_lineflux_results = {
-    u.Jy/u.sr: u.Unit('W/(m2*sr)'),
+    u.Jy/u.sr: u.Unit('W/(m2)'),
     u.Jy: u.Unit('W/m2'),
-    u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2*sr)'),
+    u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2)'),
     u.Unit('W/m2/m'): u.Unit('W/m2')
 }
 

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

@@ -200,6 +200,9 @@ def test_sb_unit_conversion(cubeviz_helper):
     viewer_1d = cubeviz_helper.app.get_viewer(
         cubeviz_helper._default_spectrum_viewer_reference_name)
 
+    uc_plg._obj.show_translator = True
+    uc_plg.flux_or_sb.selected = 'Surface Brightness'
+
     # Surface Brightness conversion
     uc_plg.flux_or_sb_unit = 'Jy / sr'
     y_display_unit = u.Unit(viewer_1d.state.y_display_unit)

jdaviz/conftest.py

@@ -214,7 +214,7 @@ def _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy, shape=(2, 2, 4), with_u
                 "TELESCOP": "JWST", "BUNIT": fluxunit.to_string(), "PIXAR_A2": 0.01}
     w = WCS(wcs_dict)
     if with_uncerts:
-        uncert = StdDevUncertainty(np.abs(np.random.normal(flux) * u.Jy))
+        uncert = StdDevUncertainty(np.abs(np.random.normal(flux) * fluxunit))
 
         return Spectrum1D(flux=flux,
                           uncertainty=uncert,
@@ -256,6 +256,12 @@ def spectrum1d_cube_custom_fluxunit():
     return _create_spectrum1d_cube_with_fluxunit
 
 
+@pytest.fixture
+def spectrum1d_cube_fluxunit_jy_per_steradian():
+    return _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy/u.sr, shape=(10, 4, 5),
+                                                 with_uncerts=True)
+
+
 @pytest.fixture
 def mos_spectrum1d(mos_spectrum2d):
     '''

jdaviz/core/helpers.py

@@ -27,8 +27,7 @@
 from jdaviz.app import Application
 from jdaviz.core.events import SnackbarMessage, ExitBatchLoadMessage
 from jdaviz.core.template_mixin import show_widget
-from jdaviz.utils import data_has_valid_wcs
-from jdaviz.app import UnitConverterWithSpectral as uc
+from jdaviz.utils import data_has_valid_wcs, flux_conversion, spectral_axis_conversion
 
 
 __all__ = ['ConfigHelper', 'ImageConfigHelper']
@@ -480,19 +479,18 @@ def _handle_display_units(self, data, use_display_units=True):
                     else:
                         # if not specified as NDUncertainty, assume stddev:
                         new_uncert = uncertainty
-                    new_uncert_converted = uc._flux_conversion(self, data,
-                                                               new_uncert.quantity.value,
-                                                               new_uncert.unit, flux_unit)
+                    new_uncert_converted = flux_conversion(data, new_uncert.quantity.value,
+                                                           new_uncert.unit, flux_unit)
                     new_uncert = StdDevUncertainty(new_uncert_converted, unit=flux_unit)
                 else:
                     new_uncert = None
 
-                new_flux = uc._flux_conversion(self, data, data.flux.value, data.flux.unit,
-                                               flux_unit) * u.Unit(flux_unit)
-                new_spec = uc._spectral_axis_conversion(self,
-                                                        data.spectral_axis.value,
-                                                        data.spectral_axis.unit,
-                                                        spectral_unit) * u.Unit(spectral_unit)
+                new_flux = flux_conversion(data, data.flux.value, data.flux.unit,
+                                           flux_unit) * u.Unit(flux_unit)
+                new_spec = (spectral_axis_conversion(data.spectral_axis.value,
+                                                     data.spectral_axis.unit,
+                                                     spectral_unit)
+                            * u.Unit(spectral_unit))
 
                 data = Spectrum1D(spectral_axis=new_spec,
                                   flux=new_flux,

jdaviz/tests/test_app.py

@@ -6,7 +6,7 @@
 from specutils import Spectrum1D
 from jdaviz import Application, Specviz
 from jdaviz.configs.default.plugins.gaussian_smooth.gaussian_smooth import GaussianSmooth
-from jdaviz.app import UnitConverterWithSpectral as uc
+from jdaviz.utils import flux_conversion
 
 
 # This applies to all viz but testing with Imviz should be enough.
@@ -227,8 +227,8 @@ def test_to_unit(cubeviz_helper):
     original_units = u.MJy / u.sr
     target_units = u.MJy
 
-    value = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
-                                original_units, target_units)
+    value = flux_conversion(data.get_object(cls=Spectrum1D),
+                            values, original_units, target_units)
 
     # will be a uniform array since not wavelength dependent
     # so test first value in array
@@ -240,8 +240,8 @@ def test_to_unit(cubeviz_helper):
     original_units = u.MJy
     target_units = u.erg / u.cm**2 / u.s / u.AA
 
-    new_values = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
-                                     original_units, target_units)
+    new_values = flux_conversion(data.get_object(cls=Spectrum1D), values,
+                                 original_units, target_units)
 
     assert np.allclose(new_values,
                        (values * original_units)
@@ -255,8 +255,8 @@ def test_to_unit(cubeviz_helper):
     original_units = u.MJy
     target_units = u.erg / u.cm**2 / u.s / u.AA
 
-    new_values = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
-                                     original_units, target_units)
+    new_values = flux_conversion(data.get_object(cls=Spectrum1D), values,
+                                 original_units, target_units)
 
     # In this case we do a regular spectral density conversion, but using the
     # first value in the spectral axis for the equivalency