convert flux cubes to per-square-pixel surface brightness cubes

CHANGES.rst

@@ -6,7 +6,7 @@ New Features
 
 - Added flux/surface brightness translation and surface brightness
   unit conversion in Cubeviz and Specviz. [#2781, #2940, #3088, #3111, #3113, #3129,
-  #3139, #3149, #3155, #3178, #3185, #3187, #3190]
+  #3139, #3149, #3155, #3178, #3185, #3187, #3190, #3156]
 
 - Plugin tray is now open by default. [#2892]
 

jdaviz/app.py

@@ -78,15 +78,27 @@
                     'Jy', 'mJy', 'uJy', 'MJy',
                     'W / (Hz m2)', 'eV / (Hz s m2)',
                     'erg / (Hz s cm2)', 'erg / (Angstrom s cm2)',
-                    'ph / (Angstrom s cm2)', 'ph / (Hz s cm2)'
+                    'ph / (Angstrom s cm2)', 'ph / (Hz s cm2)',
+                    'ct'
                 ]
                 + [
                     'Jy / sr', 'mJy / sr', 'uJy / sr', 'MJy / sr',
                     'W / (Hz sr m2)', 'eV / (Hz s sr m2)',
                     'erg / (s sr cm2)', 'erg / (Hz s sr cm2)',
                     'erg / (Angstrom s sr cm2)',
-                    'ph / (Angstrom s sr cm2)', 'ph / (Hz s sr cm2)'
-                ])
+                    'ph / (Angstrom s sr cm2)', 'ph / (Hz s sr cm2)',
+                    'ct / sr'
+                ]
+                + [
+                    'Jy / pix2', 'mJy / pix2', 'uJy / pix2', 'MJy / pix2',
+                    'W / (Hz m2 pix2)', 'eV / (Hz s m2 pix2)',
+                    'erg / (s cm2 pix2)', 'erg / (Hz s cm2 pix2)',
+                    'erg / (Angstrom s cm2 pix2)',
+                    'ph / (Angstrom s cm2 pix2)', 'ph / (Hz s cm2 pix2)',
+                    'ct / pix2'
+                ]
+
+                )
         else:  # spectral axis
             # prefer Hz over Bq and um over micron
             exclude = {'Bq', 'micron'}
@@ -1294,16 +1306,38 @@
                 sv_y_unit = sv.data()[0].flux.unit
                 if axis == 'spectral_y':
                     return sv_y_unit
-                elif axis == 'flux':
-                    if check_if_unit_is_per_solid_angle(sv_y_unit):
-                        # TODO: this will need updating once solid-angle unit can be non-steradian
-                        return sv_y_unit * u.sr
+                # since this is where we're falling back on native units (UC plugin might not exist)
+                # first check the spectrum viewer y axis for any solid angle unit (i think that it
+                # will ALWAYS be in flux, but just to be sure). If no solid angle unit is found,
+                # check the flux viewer for surface brightness units
+                sv_y_angle_unit = check_if_unit_is_per_solid_angle(sv_y_unit, return_unit=True)
+
+                # check flux viewer if none in spectral viewer
+                fv_angle_unit = None
+                if not sv_y_angle_unit:
+                    if hasattr(self._jdaviz_helper, '_default_flux_viewer_reference_name'):
+                        vname = self._jdaviz_helper._default_flux_viewer_reference_name
+                        fv = self.get_viewer(vname)
+                        fv_unit = fv.data()[0].get_object().flux.unit
+                        fv_angle_unit = check_if_unit_is_per_solid_angle(fv_unit,
+                                                                         return_unit=True)
+                    else:
+                        # mosviz, not sure what to do here but can't access flux
+                        # viewer the same way. once we force the UC plugin to
+                        # exist this won't matter anyway because units can be
+                        # acessed from the plugin directly. assume u.sr for now
+                        fv_angle_unit = u.sr
+
+                solid_angle_unit = sv_y_angle_unit or fv_angle_unit
+
+                if axis == 'flux':
+                    if sv_y_angle_unit:
+                        return sv_y_unit * solid_angle_unit
                     return sv_y_unit
-                else:
-                    # surface_brightness
-                    if check_if_unit_is_per_solid_angle(sv_y_unit):
+                elif axis == 'sb':
+                    if sv_y_angle_unit:
                         return sv_y_unit
-                    return sv_y_unit / u.sr
+                    return sv_y_unit / solid_angle_unit
             elif axis == 'temporal':
                 # No unit for ramp's time (group/resultant) axis:
                 return None

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

@@ -20,7 +20,6 @@
                                         SpectralContinuumMixin,
                                         skip_if_no_updates_since_last_active,
                                         with_spinner)
-from jdaviz.core.validunits import check_if_unit_is_per_solid_angle
 from jdaviz.core.user_api import PluginUserApi
 
 __all__ = ['MomentMap']
@@ -111,8 +110,10 @@ def __init__(self, *args, **kwargs):
         self.output_unit = SelectPluginComponent(self,
                                                  items='output_unit_items',
                                                  selected='output_unit_selected',
-                                                 manual_options=['Flux', 'Spectral Unit',
-                                                                 'Velocity', 'Velocity^N'])
+                                                 manual_options=['Surface Brightness',
+                                                                 'Spectral Unit',
+                                                                 'Velocity',
+                                                                 'Velocity^N'])
 
         self.dataset.add_filter('is_cube')
         self.add_results.viewer.filters = ['is_image_viewer']
@@ -174,16 +175,11 @@ def _set_data_units(self, event={}):
             self.output_unit_selected = moment_unit_options[unit_options_index][0]
         self.send_state("output_unit_selected")
 
-        # either 'Flux' or 'Surface Brightness'
-        orig_spectral_y_type = self.output_unit_items[0]['label']
-
-        unit_dict = {orig_spectral_y_type: "",
+        unit_dict = {"Surface Brightness": "",
                      "Spectral Unit": "",
                      "Velocity": "km/s",
                      "Velocity^N": f"km{self.n_moment}/s{self.n_moment}"}
 
-        sb_or_flux_label = None
-
         if self.dataset_selected != "":
             # Spectral axis is first in this list
             data = self.app.data_collection[self.dataset_selected]
@@ -196,37 +192,11 @@ def _set_data_units(self, event={}):
                 sunit = ""
             self.dataset_spectral_unit = sunit
             unit_dict["Spectral Unit"] = sunit
-
-            # get flux/SB units
-            if self.spectrum_viewer and hasattr(self.spectrum_viewer.state, 'y_display_unit'):
-                if self.spectrum_viewer.state.y_display_unit is not None:
-                    unit_dict[orig_spectral_y_type] = self.spectrum_viewer.state.y_display_unit
-                else:
-                    # spectrum_viewer.state will only have x/y_display_unit if unit conversion has
-                    # been done if not, get default flux units which should be the units displayed
-                    unit_dict[orig_spectral_y_type] = data.get_component('flux').units
-            else:
-                # spectrum_viewer.state will only have x/y_display_unit if unit conversion has
-                # been done if not, get default flux units which should be the units displayed
-                unit_dict[orig_spectral_y_type] = data.get_component('flux').units
-
-            # figure out if label should say 'Flux' or 'Surface Brightness'
-            sb_or_flux_label = "Flux"
-            is_unit_solid_angle = check_if_unit_is_per_solid_angle(unit_dict[orig_spectral_y_type])
-            if is_unit_solid_angle is True:
-                sb_or_flux_label = "Surface Brightness"
+            unit_dict["Surface Brightness"] = self.app._get_display_unit('sb')
 
         # Update units in selection item dictionary
         for item in self.output_unit_items:
             item["unit_str"] = unit_dict[item["label"]]
-            if item["label"] in ["Flux", "Surface Brightness"] and sb_or_flux_label:
-                # change unit label to reflect if unit is flux or SB
-                item["label"] = sb_or_flux_label
-
-        if self.dataset_selected != "":
-            # output_unit_selected might not match (potentially) changed flux/SB label
-            if self.output_unit_selected in ['Flux', 'Surface Brightness']:
-                self.output_unit_selected = sb_or_flux_label
 
         # Filter what we want based on n_moment
         if self.n_moment == 0:
@@ -343,18 +313,18 @@ def calculate_moment(self, add_data=True):
         # convert units for moment 0, which is the only currently supported
         # moment for using converted units.
         if n_moment == 0:
-            # get flux/SB units
-            if self.spectrum_viewer and hasattr(self.spectrum_viewer.state, 'y_display_unit'):
-                if self.spectrum_viewer.state.y_display_unit is not None:
-                    flux_sb_unit = self.spectrum_viewer.state.y_display_unit
-                else:
-                    flux_sb_unit = data.get_component('flux').units
-            else:
-                flux_sb_unit = data.get_component('flux').units
+            # get display units for moment 0 based on unit conversion plugin selection
+            # the 0th item of this dictionary is the surface brightness unit, kept
+            # up to date with choices from the UC plugin
+            moment_0_display_unit = self.output_unit_items[0]['unit_str']
+
+            # convert unit string to Unit so moment map data can be converted
+            # `display_unit` is the data unit (surface brighntess) translated
+            # to the choice of selected flux and angle unit from the UC plugin.
+            display_unit = u.Unit(moment_0_display_unit)
+
+            moment_new_unit = display_unit * self.app._get_display_unit('spectral')  # noqa: E501
 
-            # convert unit string to u.Unit so moment map data can be converted
-            spectral_y_display_unit = u.Unit(flux_sb_unit)
-            moment_new_unit = spectral_y_display_unit * self.spectrum_viewer.state.x_display_unit  # noqa: E501
             self.moment = self.moment.to(moment_new_unit)
 
         # Reattach the WCS so we can load the result

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

@@ -12,10 +12,21 @@
 from numpy.testing import assert_allclose
 
 
-def test_user_api(cubeviz_helper, spectrum1d_cube):
+@pytest.mark.parametrize("cube_type", ["Surface Brightness", "Flux"])
+def test_user_api(cubeviz_helper, spectrum1d_cube, spectrum1d_cube_sb_unit, cube_type):
+
+    # test is parameterize to test a cube that is in Jy / sr (Surface Brightness)
+    # as well as Jy (Flux), to test that flux cubes, which are converted in the
+    # parser to flux / pix^2 surface brightness cubes, both work correctly.
+
+    if cube_type == 'Surface Brightness':
+        cube = spectrum1d_cube_sb_unit
+    elif cube_type == 'Flux':
+        cube = spectrum1d_cube
+
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="No observer defined on WCS.*")
-        cubeviz_helper.load_data(spectrum1d_cube, data_label='test')
+        cubeviz_helper.load_data(cube, data_label='test')
 
     mm = cubeviz_helper.plugins['Moment Maps']
     assert not mm._obj.continuum_marks['center'].visible
@@ -43,20 +54,33 @@ def test_user_api(cubeviz_helper, spectrum1d_cube):
         mm.n_moment = 1
         with pytest.raises(ValueError, match="not one of"):
             mm._obj.output_unit_selected = "Bad Unit"
-        mm._obj.output_unit_selected = "Flux"
+        mm._obj.output_unit_selected = "Surface Brightness"
         with pytest.raises(ValueError, match="units must be in"):
             mm.calculate_moment()
 
 
-def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path):
+@pytest.mark.parametrize("cube_type", ["Surface Brightness", "Flux"])
+def test_moment_calculation(cubeviz_helper, spectrum1d_cube,
+                            spectrum1d_cube_sb_unit, cube_type, tmp_path):
 
     moment_unit = "Jy m"
     moment_value_str = "+6.40166e-10"
 
+    if cube_type == 'Surface Brightness':
+        moment_unit += " / sr"
+        cube = spectrum1d_cube_sb_unit
+        cube_unit = cube.unit.to_string()
+
+    elif cube_type == 'Flux':
+        moment_unit += " / pix2"
+        cube = spectrum1d_cube
+        cube_unit = cube.unit.to_string() + " / pix2"  # cube in Jy will become cube in Jy / pix2
+
     dc = cubeviz_helper.app.data_collection
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="No observer defined on WCS.*")
-        cubeviz_helper.load_data(spectrum1d_cube, data_label='test')
+        cubeviz_helper.load_data(cube, data_label='test')
+
     flux_viewer = cubeviz_helper.app.get_viewer(cubeviz_helper._default_flux_viewer_reference_name)
 
     # Since we are not really displaying, need this to trigger GUI stuff.
@@ -68,7 +92,9 @@ def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path):
 
     mm.n_moment = 0  # Collapsed sum, will get back 2D spatial image
     assert mm._obj.results_label == 'moment 0'
-    assert mm.output_unit == "Flux"
+
+    # result is always a SB unit - if flux cube loaded, per pix2
+    assert mm.output_unit == 'Surface Brightness'
 
     mm._obj.add_results.viewer.selected = cubeviz_helper._default_uncert_viewer_reference_name
     mm._obj.vue_calculate_moment()
@@ -88,13 +114,13 @@ def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path):
     assert mm._obj.results_label == 'moment 0'
     assert mm._obj.results_label_overwrite is True
 
-    # Make sure coordinate display works
+    # Make sure coordinate display works in flux viewer (loaded data, not the moment map)
     label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info']
     label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove',
                                                       'domain': {'x': 0, 'y': 0}})
     assert flux_viewer.state.slices == (0, 0, 1)
     # Slice 0 has 8 pixels, this is Slice 1
-    assert label_mouseover.as_text() == ("Pixel x=00.0 y=00.0 Value +8.00000e+00 Jy",
+    assert label_mouseover.as_text() == (f"Pixel x=00.0 y=00.0 Value +8.00000e+00 {cube_unit}",
                                          "World 13h39m59.9731s +27d00m00.3600s (ICRS)",
                                          "204.9998877673 27.0001000000 (deg)")
 
@@ -290,7 +316,7 @@ def test_momentmap_nirspec_prism(cubeviz_helper, tmp_path):
                     (sky_cube.ra.deg, sky_cube.dec.deg))
 
 
-def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_fluxunit):
+def test_correct_output_spectral_y_units(cubeviz_helper, spectrum1d_cube_custom_fluxunit):
 
     moment_unit = "Jy m / sr"