Skip to content

Commit

Permalink
Fix display units of moments in velocity (#2665)
Browse files Browse the repository at this point in the history
* 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>
  • Loading branch information
rosteen and pllim authored Jan 18, 2024
1 parent 5c69d28 commit 9140858
Showing 3 changed files with 23 additions and 18 deletions.
2 changes: 1 addition & 1 deletion CHANGES.rst
Original file line number Diff line number Diff line change
@@ -13,7 +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]

28 changes: 16 additions & 12 deletions jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
Original file line number Diff line number Diff line change
@@ -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"
Original file line number Diff line number Diff line change
@@ -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):

0 comments on commit 9140858

Please sign in to comment.