Fix unit handling in continuum and moment calculations (spacetelescop…

…e#3216)

* Fix unit handling in continuum and moment calculations

* Add test that fails on main

* Bump specutils pin

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, #3156, #3200, #3192, #3206, #3211]
+  #3139, #3149, #3155, #3178, #3185, #3187, #3190, #3156, #3200, #3192, #3206, #3211, #3216]
 
 - Plugin tray is now open by default. [#2892]
 

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

@@ -271,7 +271,13 @@ def calculate_moment(self, add_data=True):
 
         # slice out desired region
         # TODO: should we add a warning for a composite spectral subset?
-        spec_reg = self.spectral_subset.selected_obj
+        if self.spectral_subset.selected == "Entire Spectrum":
+            spec_reg = None
+        else:
+            spec_reg = self.app.get_subsets(self.spectral_subset.selected,
+                                            simplify_spectral=True,
+                                            use_display_units=True)
+        # We need to convert the spectral region to the display units
 
         if spec_reg is None:
             slab = cube
@@ -303,7 +309,12 @@ def calculate_moment(self, add_data=True):
             ref_wavelength = self.reference_wavelength * u.Unit(self.dataset_spectral_unit)
             slab_sa = slab.spectral_axis.to("km/s", doppler_convention="relativistic",
                                             doppler_rest=ref_wavelength)
-            slab = Spectrum1D(slab.flux, slab_sa)
+            slab = Spectrum1D(slab.flux, slab_sa, uncertainty=slab.uncertainty)
+        # Otherwise convert spectral axis to display units, have to do frequency <-> wavelength
+        # before calculating
+        else:
+            slab_sa = slab.spectral_axis.to(self.app._get_display_unit('spectral'))
+            slab = Spectrum1D(slab.flux, slab_sa, uncertainty=slab.uncertainty)
 
         # Finally actually calculate the moment
         self.moment = analysis.moment(slab, order=n_moment).T

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

@@ -7,6 +7,7 @@
 from astropy import units as u
 from astropy.io import fits
 from astropy.nddata import CCDData
+from astropy.tests.helper import assert_quantity_allclose
 from astropy.wcs import WCS
 from numpy.testing import assert_allclose
 from specutils import SpectralRegion
@@ -244,9 +245,12 @@ def test_moment_frequency_unit_conversion(cubeviz_helper, spectrum1d_cube_larger
     mm.n_moment = 1
     mm.output_unit = 'Spectral Unit'
     moment_1_data = mm.calculate_moment()
+    mm.n_moment = 0
+    moment_0_data = mm.calculate_moment()
 
     # Check to make sure there are no nans
     assert len(np.where(moment_1_data.data > 0)[0]) == 8
+    assert_quantity_allclose(moment_0_data, -2.9607526e+09*u.Unit("Hz Jy / pix2"))
 
 
 def test_write_momentmap(cubeviz_helper, spectrum1d_cube, tmp_path):

jdaviz/core/template_mixin.py

@@ -2970,7 +2970,8 @@ def _get_continuum(self, dataset, spectral_subset, update_marks=False, per_pixel
         if per_pixel:
             if self.app.config != 'cubeviz':
                 raise ValueError("per-pixel only supported for cubeviz")
-            full_spectrum = self.app._jdaviz_helper.get_data(self.dataset.selected)
+            full_spectrum = self.app._jdaviz_helper.get_data(self.dataset.selected,
+                                                             use_display_units=True)
         else:
             full_spectrum = dataset.get_selected_spectrum(use_display_units=True)
 
@@ -2994,7 +2995,8 @@ def _get_continuum(self, dataset, spectral_subset, update_marks=False, per_pixel
             spectrum = full_spectrum
         else:
             sr = self.app.get_subsets(spectral_subset.selected,
-                                      simplify_spectral=True)
+                                      simplify_spectral=True,
+                                      use_display_units=True)
             spectrum = extract_region(full_spectrum, sr, return_single_spectrum=True)
             sr_lower = np.nanmin(spectrum.spectral_axis[spectrum.spectral_axis >= sr.lower])  # noqa
             sr_upper = np.nanmax(spectrum.spectral_axis[spectrum.spectral_axis <= sr.upper])  # noqa

pyproject.toml

@@ -23,7 +23,7 @@ dependencies = [
     "ipywidgets>=8.0.6",
     "solara>=1.40.0",
     "pyyaml>=5.4.1",
-    "specutils>=1.16",
+    "specutils>=1.18",
     "specreduce>=1.4.1",
     "photutils>=1.4",
     "glue-astronomy>=0.10",