convert flux cubes to per-square-pixel surface brightness cubes (#3156)

* checkpoint

* clean up

* code style

* remove debugging

* fixed equivs, more tests passing?

* fix untranslatable units for per pixels

* fixed mm

* .

* fix mouseover for sr

* fixed batch ap phot

* fix lingering issues

* change log

* code style

* doc

* checkpoint

* removing example from docstring because it keeps causing issues

* line analysis fixes

* code style

* fix docs

* test

* .

* reverting wrong change selected in rebase

* more missing code from main

* review comment

* small review changes

* code style

* fix formatting of unit after appending angle

* review comments, skip test

* review comments, counts support

* code style

* added test, cleanup

* skip test in viewers

* skip test in parsers

* remove var def to avoid line overflow

* remove new traitlet in ap_phot

* code style
;

* added back in old test, updated docstrings

* review comments

* pllim review from 3192

* update for imviz aperphot

* final review comments from kyle

* Fix NaN handling in cube fitting and initial fixes for unit conversion/model fitting interaction

* Remove debugging prints, add comment for context

* Codestyle, changelog

* Reestimate parameters when cube fitting is toggled

* Only reestimate if spectral y type isn't SB

* Codestyle

* Changelog

* Fix initializing linear component for cube fit

* Handle linear component estimation for cube case

* Only reshape here in 3D case

* Skip tests that need #3156

* Respect selected display units when initializing model components

* Use app._get_display_unit instead of relying on Unit Conversion

* Add equivalency here

* Only reestimate here if sb unit != spectral_y unit

* Don't automatically reestimate when toggling cube fit, make the user do it

* Remove print

* Check for warning in test after cube toggle

Fix test

Codestyle

* Add test for cube fitting after flux unit change

* Add a to-do about a test for unit conversion with equivalency

* Fix failing test

* review comments

* .

* model fitting tests

* Back to parallel processing post-debugging

* review comments

* code style

* move import

---------

Co-authored-by: Ricky O'Steen <[email protected]>

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 @@ def equivalent_units(self, data, cid, units):
  '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 @@ def _get_display_unit(self, axis):
  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"