Skip to content

Commit

Permalink
Add test that triggers handle_display_units bug (spacetelescope#2910)
Browse files Browse the repository at this point in the history 
* Add test that triggers handle_display_units bug

* Pass test with hardcoded values

* Pass style checks

* Update handle_display_units

* Fix style

* Fix import and style issues

* Rename variable

* Try to handle 2D objects using Spectrum1D

* Remove commented code

* Fix moment map with translator on, leave comments for continuum

* Fix code style

* Remove wcs from NDDataArray in to_unit

* Revert review changes

* Comment out test until unit conversion translator is active

* Pass wcs when converting to spectrum1d from nddataarray

* Do not pass wcs

* Add change log
  • Loading branch information
@javerbukh
javerbukh authored Jun 18, 2024
1 parent e0b1ec3 commit a1d9682
Show file tree
Hide file tree
Showing 8 changed files with 139 additions and 76 deletions.
2 changes: 2 additions & 0 deletions CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,8 @@ Cubeviz
- Update the scale factor used to convert a spectrum between surface brightness and flux
to use wavelength-dependent aperture area instead of the cone slice scale factor. [#2860]

- Handle display units when doing flux / surface brightness conversions. [#2910]

Imviz
^^^^^

Expand Down
131 changes: 67 additions & 64 deletions jdaviz/app.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import warnings
import ipyvue
from astropy import units as u
from astropy.nddata import NDData
from astropy.nddata import NDData, NDDataArray
from astropy.io import fits
from astropy.time import Time
from astropy.units import Quantity
Expand Down Expand Up @@ -53,7 +53,7 @@
from jdaviz.utils import (SnackbarQueue, alpha_index, data_has_valid_wcs, layer_is_table_data,
MultiMaskSubsetState, _wcs_only_label)

__all__ = ['Application', 'ALL_JDAVIZ_CONFIGS']
__all__ = ['Application', 'ALL_JDAVIZ_CONFIGS', 'UnitConverterWithSpectral']

SplitPanes()
GoldenLayout()
Expand Down Expand Up @@ -109,72 +109,75 @@ def to_unit(self, data, cid, values, original_units, target_units):
spec = data.get_object(cls=Spectrum1D)

except RuntimeError:
eqv = []
data = data.get_object(cls=NDDataArray)
spec = Spectrum1D(flux=data.data * u.Unit(original_units))
return self._flux_conversion(spec, values, original_units, target_units)
else: # spectral axis
return self._spectral_axis_conversion(values, original_units, target_units)

def _flux_conversion(self, spec, values, original_units, target_units):
# If there are only two values, this is likely the limits being converted, so then
# in case we need to use the spectral density equivalency, we need to provide only
# to spectral axis values. If there is only one value
if not np.isscalar(values) and len(values) == 2:
spectral_values = spec.spectral_axis[0]
else:
spectral_values = spec.spectral_axis

# Ensure a spectrum passed through Spectral Extraction plugin
if '_pixel_scale_factor' in spec.meta and len(values) != 2:
# Data item in data collection does not update from conversion/translation.
# App wide original data units are used for conversion, original and
# target_units dictate the conversion to take place.

if (u.sr in u.Unit(original_units).bases) and \
(u.sr not in u.Unit(target_units).bases):
# Surface Brightness -> Flux
eqv = [(u.MJy / u.sr,
u.MJy,
lambda x: (x * np.array(spec.meta.get('_pixel_scale_factor', 1))),
lambda x: x)]
elif (u.sr not in u.Unit(original_units).bases) and \
(u.sr in u.Unit(target_units).bases):
# Flux -> Surface Brightness
eqv = [(u.MJy,
u.MJy / u.sr,
lambda x: (x / np.array(spec.meta.get('_pixel_scale_factor', 1))),
lambda x: x)]
else:
eqv = []

# If there are only two values, this is likely the limits being converted, so then
# in case we need to use the spectral density equivalency, we need to provide only
# to spectral axis values. If there is only one value
if not np.isscalar(values) and len(values) == 2:
spectral_values = spec.spectral_axis[0]
else:
spectral_values = spec.spectral_axis

# Ensure a spectrum passed through Spectral Extraction plugin
if '_pixel_scale_factor' in spec.meta and len(values) != 2:

# Data item in data collection does not update from conversion/translation.
# App wide original data units are used for conversion, original and
# target_units dictate the conversion to take place.

if (u.sr in u.Unit(original_units).bases) and \
(u.sr not in u.Unit(target_units).bases):
# Surface Brightness -> Flux
eqv = [(u.MJy / u.sr,
u.MJy,
lambda x: (x * np.array(spec.meta['_pixel_scale_factor'])),
lambda x: x)]
elif (u.sr not in u.Unit(original_units).bases) and \
(u.sr in u.Unit(target_units).bases):
# Flux -> Surface Brightness
eqv = [(u.MJy,
u.MJy / u.sr,
lambda x: (x / np.array(spec.meta['_pixel_scale_factor'])),
lambda x: x)]
else:
eqv = u.spectral_density(spectral_values)

elif len(values) == 2:
# Need this for setting the y-limits
eqv = u.spectral_density(spectral_values)

if '_pixel_scale_factor' in spec.meta:
# get min and max scale factors, to use with min and max of spec for
# y-limits. Make sure they are Quantities (can be numpy.float64).
pixel_scale_min = (Quantity(min(spec.meta['_pixel_scale_factor']))).value
pixel_scale_max = (Quantity(max(spec.meta['_pixel_scale_factor']))).value
min_max = [pixel_scale_min, pixel_scale_max]

if (u.sr in u.Unit(original_units).bases) and \
(u.sr not in u.Unit(target_units).bases):
eqv += [(u.MJy,
u.MJy / u.sr,
lambda x: x * np.array(min_max),
lambda x: x)]
elif (u.sr not in u.Unit(original_units).bases) and \
(u.sr in u.Unit(target_units).bases):
eqv += [(u.MJy / u.sr,
u.MJy,
lambda x: x / np.array(min_max),
lambda x: x)]
eqv = u.spectral_density(spectral_values)

elif len(values) == 2:
# Need this for setting the y-limits
eqv = u.spectral_density(spectral_values)

if '_pixel_scale_factor' in spec.meta:
# get min and max scale factors, to use with min and max of spec for
# y-limits. Make sure they are Quantities (can be numpy.float64).
pixel_scale_min = (Quantity(min(spec.meta.get('_pixel_scale_factor', 1)))).value
pixel_scale_max = (Quantity(max(spec.meta.get('_pixel_scale_factor', 1)))).value
min_max = [pixel_scale_min, pixel_scale_max]

if (u.sr in u.Unit(original_units).bases) and \
(u.sr not in u.Unit(target_units).bases):
eqv += [(u.MJy,
u.MJy / u.sr,
lambda x: x * np.array(min_max),
lambda x: x)]
elif (u.sr not in u.Unit(original_units).bases) and \
(u.sr in u.Unit(target_units).bases):
eqv += [(u.MJy / u.sr,
u.MJy,
lambda x: x / np.array(min_max),
lambda x: x)]

else:
eqv = u.spectral_density(spectral_values)
else:
eqv = u.spectral_density(spectral_values)

else: # spectral axis
eqv = u.spectral() + u.pixel_scale(1*u.pix)
return (values * u.Unit(original_units)).to_value(u.Unit(target_units), equivalencies=eqv)

def _spectral_axis_conversion(self, values, original_units, target_units):
eqv = u.spectral() + u.pixel_scale(1*u.pix)
return (values * u.Unit(original_units)).to_value(u.Unit(target_units), equivalencies=eqv)


Expand Down
15 changes: 14 additions & 1 deletion jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@
with_spinner)
from jdaviz.core.validunits import check_if_unit_is_per_solid_angle
from jdaviz.core.user_api import PluginUserApi
from jdaviz.app import UnitConverterWithSpectral as uc


__all__ = ['MomentMap']

Expand Down Expand Up @@ -356,7 +358,18 @@ def calculate_moment(self, add_data=True):
moment_new_unit = flux_or_sb_display_unit
else:
moment_new_unit = flux_or_sb_display_unit * self.spectrum_viewer.state.x_display_unit # noqa: E501
self.moment = self.moment.to(moment_new_unit)

# Create a temporary Spectrum1D object with ability to convert from surface brightness
# to flux
temp_spec = Spectrum1D(flux=self.moment)
flux_values = np.sum(np.ones_like(temp_spec.flux.value), axis=(0, 1))
pix_scale = self.dataset.selected_dc_item.meta.get('PIXAR_SR', 1.0)
pix_scale_factor = (flux_values * pix_scale)
temp_spec.meta['_pixel_scale_factor'] = pix_scale_factor
converted_spec = uc._flux_conversion(self, temp_spec, self.moment.value,
self.moment.unit,
moment_new_unit) * moment_new_unit
self.moment = converted_spec

# Reattach the WCS so we can load the result
self.moment = CCDData(self.moment, wcs=data_wcs)
Expand Down
19 changes: 19 additions & 0 deletions jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,3 +334,22 @@ def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_
# and that calculated moment has the correct units
mm.calculate_moment()
assert mm.moment.unit == moment_unit

uc._obj.show_translator = True
uc.flux_or_sb.selected = 'Flux'
mm._set_data_units()

# and make sure this change is propogated
output_unit_moment_0 = mm.output_unit_items[0]
assert output_unit_moment_0['label'] == 'Flux'
assert output_unit_moment_0['unit_str'] == 'Jy'

# TODO: Failing because of dev version of upstream dependency, figure
# out which one
# assert mm.calculate_moment()

# TODO: This test should pass once continuum subtraction works with
# flux to surface brightness conversion
# mm.continuum.selected = 'Surrounding'
#
# assert mm.calculate_moment()
3 changes: 3 additions & 0 deletions jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,3 +232,6 @@ def test_sb_unit_conversion(cubeviz_helper):
y_display_unit = u.Unit(viewer_1d.state.y_display_unit)

assert y_display_unit == u.MJy

la = cubeviz_helper.plugins['Line Analysis']._obj
assert la.dataset.get_selected_spectrum(use_display_units=True)
24 changes: 17 additions & 7 deletions jdaviz/core/helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@
from jdaviz.core.events import SnackbarMessage, ExitBatchLoadMessage
from jdaviz.core.template_mixin import show_widget
from jdaviz.utils import data_has_valid_wcs
from jdaviz.app import UnitConverterWithSpectral as uc


__all__ = ['ConfigHelper', 'ImageConfigHelper']

Expand Down Expand Up @@ -474,18 +476,26 @@ def _handle_display_units(self, data, use_display_units=True):
if hasattr(uncertainty, 'represent_as'):
new_uncert = uncertainty.represent_as(
StdDevUncertainty
).quantity.to(flux_unit)
)
else:
# if not specified as NDUncertainty, assume stddev:
new_uncert = uncertainty.quantity.to(flux_unit)
new_uncert = StdDevUncertainty(new_uncert, unit=flux_unit)
new_uncert = uncertainty
new_uncert_converted = uc._flux_conversion(self, data,
new_uncert.quantity.value,
new_uncert.unit, flux_unit)
new_uncert = StdDevUncertainty(new_uncert_converted, unit=flux_unit)
else:
new_uncert = None

data = Spectrum1D(spectral_axis=data.spectral_axis.to(spectral_unit,
u.spectral()),
flux=data.flux.to(flux_unit,
u.spectral_density(data.spectral_axis)),
new_flux = uc._flux_conversion(self, data, data.flux.value, data.flux.unit,
flux_unit) * u.Unit(flux_unit)
new_spec = uc._spectral_axis_conversion(self,
data.spectral_axis.value,
data.spectral_axis.unit,
spectral_unit) * u.Unit(spectral_unit)

data = Spectrum1D(spectral_axis=new_spec,
flux=new_flux,
uncertainty=new_uncert,
mask=data.mask)
else: # pragma: nocover
Expand Down
11 changes: 11 additions & 0 deletions jdaviz/core/template_mixin.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2946,6 +2946,17 @@ def _get_continuum(self, dataset, spectral_subset, update_marks=False, per_pixel
raise ValueError("per-pixel only supported for cubeviz")
full_spectrum = self.app._jdaviz_helper.get_data(self.dataset.selected,
use_display_units=True)
# TODO: Something like the following code may be needed to get continuum
# with display units working
# temp_spec = self.app._jdaviz_helper.get_data(self.dataset.selected,
# use_display_units=True)
# flux_values = np.sum(np.ones_like(temp_spec.flux.value), axis=(0, 1))
# pix_scale = self.dataset.selected_dc_item.meta.get('PIXAR_SR', 1.0)
# pix_scale_factor = (flux_values * pix_scale)
# temp_spec.meta['_pixel_scale_factor'] = pix_scale_factor
# full_spectrum = self._specviz_helper._handle_display_units(temp_spec,
# use_display_units=True)

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

Expand Down
10 changes: 6 additions & 4 deletions jdaviz/tests/test_app.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,13 +222,13 @@ def test_to_unit(cubeviz_helper):

extract_plg.extract()

cid = cubeviz_helper.app.data_collection[0].data.find_component_id('flux')
data = cubeviz_helper.app.data_collection[-1].data
values = [1]
original_units = u.MJy / u.sr
target_units = u.MJy

value = uc.to_unit(cubeviz_helper, data, cid, values, original_units, target_units)
value = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
original_units, target_units)

# will be a uniform array since not wavelength dependent
# so test first value in array
Expand All @@ -240,7 +240,8 @@ def test_to_unit(cubeviz_helper):
original_units = u.MJy
target_units = u.erg / u.cm**2 / u.s / u.AA

new_values = uc.to_unit(cubeviz_helper, data, cid, values, original_units, target_units)
new_values = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
original_units, target_units)

assert np.allclose(new_values,
(values * original_units)
Expand All @@ -254,7 +255,8 @@ def test_to_unit(cubeviz_helper):
original_units = u.MJy
target_units = u.erg / u.cm**2 / u.s / u.AA

new_values = uc.to_unit(cubeviz_helper, data, cid, values, original_units, target_units)
new_values = uc._flux_conversion(uc, data.get_object(cls=Spectrum1D), values,
original_units, target_units)

# In this case we do a regular spectral density conversion, but using the
# first value in the spectral axis for the equivalency
Expand Down

0 comments on commit a1d9682

Please sign in to comment.