Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix unit handling in continuum and moment calculations #3216

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down
15 changes: 13 additions & 2 deletions jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down
6 changes: 4 additions & 2 deletions jdaviz/core/template_mixin.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Comment on lines -2973 to +2974
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like @cshanahan1 had just removed this. Was there a reason for that or did it just not seem to affect that use-case (see https://github.com/spacetelescope/jdaviz/pull/3211/files#r1787615549)

else:
full_spectrum = dataset.get_selected_spectrum(use_display_units=True)

Expand All @@ -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
Expand Down