Fix display units of moments in velocity (spacetelescope#2665)

* Fixing units in result, trying to debug negative second moments

* Remove unneeded import

* Update test

* Add changelog

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

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

* Fix changelog entry

---------

Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com>

CHANGES.rst

@@ -13,8 +13,7 @@ New Features
 
 Cubeviz
 ^^^^^^^
-
-- Calculated moments can now be output in velocity units. [#2584, #2588]
+- Calculated moments can now be output in velocity units. [#2584, #2588, #2665]
 
 - Added functionality to Collapse and Spectral Extraction plugins to save results to FITS file. [#2586]
 

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

@@ -2,12 +2,11 @@
 from pathlib import Path
 
 from astropy import units as u
-from astropy.constants import c
 from astropy.nddata import CCDData
 import numpy as np
 
 from traitlets import Bool, List, Unicode, observe
-from specutils import manipulation, analysis
+from specutils import manipulation, analysis, Spectrum1D
 
 from jdaviz.core.custom_traitlets import IntHandleEmpty, FloatHandleEmpty
 from jdaviz.core.events import SnackbarMessage
@@ -243,20 +242,25 @@ def calculate_moment(self, add_data=True):
         data_wcs = getattr(cube.wcs, 'celestial', None)
         if data_wcs:
             data_wcs = data_wcs.swapaxes(0, 1)  # We also transpose WCS to match.
-        self.moment = CCDData(analysis.moment(slab, order=n_moment).T, wcs=data_wcs)
-        if n_moment > 0 and self.output_unit_selected.lower()[0:8] == "velocity":
+
+        # Convert spectral axis to velocity units if desired output is in velocity
+        if n_moment > 0 and self.output_unit_selected.lower().startswith("velocity"):
             # Catch this if called from API
             if not self.reference_wavelength > 0.0:
                 raise ValueError("reference_wavelength must be set for output in velocity units.")
-            power_unit = f"{self.dataset_spectral_unit}{self.n_moment}"
-            self.moment = np.power(self.moment.convert_unit_to(power_unit), 1/self.n_moment)
-            self.moment = self.moment << u.Unit(self.dataset_spectral_unit)
+
             ref_wavelength = self.reference_wavelength * u.Unit(self.dataset_spectral_unit)
-            relative_wavelength = (self.moment-ref_wavelength)/ref_wavelength
-            in_velocity = c*relative_wavelength
-            if self.output_unit_selected.lower() == "velocity^n":
-                in_velocity = np.power(in_velocity, self.n_moment)
-            self.moment = CCDData(in_velocity, wcs=data_wcs)
+            slab_sa = slab.spectral_axis.to("km/s", doppler_convention="relativistic",
+                                            doppler_rest=ref_wavelength)
+            slab = Spectrum1D(slab.flux, slab_sa)
+
+        # Finally actually calculate the moment
+        self.moment = analysis.moment(slab, order=n_moment).T
+        # If n>1 and velocity is desired, need to take nth root of result
+        if n_moment > 0 and self.output_unit_selected.lower() == "velocity":
+            self.moment = np.power(self.moment, 1/self.n_moment)
+        # Reattach the WCS so we can load the result
+        self.moment = CCDData(self.moment, wcs=data_wcs)
 
         fname_label = self.dataset_selected.replace("[", "_").replace("]", "")
         self.filename = f"moment{n_moment}_{fname_label}.fits"

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

@@ -165,7 +165,7 @@ def test_moment_velocity_calculation(cubeviz_helper, spectrum1d_cube, tmpdir):
     label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info']
     label_mouseover._viewer_mouse_event(uncert_viewer, {'event': 'mousemove',
                                                         'domain': {'x': 0, 'y': 0}})
-    assert label_mouseover.as_text() == ("Pixel x=00.0 y=00.0 Value -4.14382e+05 m / s",
+    assert label_mouseover.as_text() == ("Pixel x=00.0 y=00.0 Value -4.14668e+02 km / s",
                                          "World 13h39m59.9731s +27d00m00.3600s (ICRS)",
                                          "204.9998877673 27.0001000000 (deg)")
 
@@ -174,10 +174,11 @@ def test_moment_velocity_calculation(cubeviz_helper, spectrum1d_cube, tmpdir):
     mm.calculate_moment()
 
     label_mouseover._viewer_mouse_event(uncert_viewer, {'event': 'mousemove',
-                                                        'domain': {'x': 0, 'y': 0}})
-    assert label_mouseover.as_text() == ("Pixel x=00.0 y=00.0 Value -2.99792e+08 m / s",
-                                         "World 13h39m59.9731s +27d00m00.3600s (ICRS)",
-                                         "204.9998877673 27.0001000000 (deg)")
+                                                        'domain': {'x': 1, 'y': 1}})
+
+    assert label_mouseover.as_text() == ("Pixel x=01.0 y=01.0 Value +2.32415e+01 km / s",
+                                         "World 13h39m59.9461s +27d00m00.7200s (ICRS)",
+                                         "204.9997755344 27.0001999998 (deg)")
 
 
 def test_write_momentmap(cubeviz_helper, spectrum1d_cube, tmp_path):