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

convert limits on spectrum y-unit change #3335

Merged
merged 15 commits into from
Dec 20, 2024
Merged
2 changes: 2 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ New Features

* Improve performance while importing multiple regions. [#3321]

* Changing flux/SB display units no longer resets viewer zoom levels. [#3335]

Cubeviz
^^^^^^^

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -356,3 +356,30 @@ def test_cubeviz_flux_sb_translation_counts(cubeviz_helper, angle_unit):
f"Pixel x=00010.0 y=00008.0 Value +1.00000e+00 ct / {angle_str}",
"World 13h39m59.7037s +27d00m03.2400s (ICRS)",
"204.9987654313 27.0008999946 (deg)")


@pytest.mark.parametrize("start_unit, end_unit, end_spectral_y_type, expected_limits",
[('Jy', 'MJy', 'Flux', (5e-07, 6e-07, 1e-4, 1.05e-4)),
('MJy', 'ph / (Angstrom s cm2)', 'Surface Brightness', (5e-07, 6e-07, 25153169.66070254, 31692993.772485193)) # noqa
])
def test_limits_on_unit_change(cubeviz_helper, start_unit, end_unit,
end_spectral_y_type, expected_limits):
"""
Test that the limits are reset when changing units
"""

w, wcs_dict = cubeviz_wcs_dict()
flux = np.zeros((30, 20, 3001), dtype=np.float32)
flux[5:15, 1:11, :] = 1
cube = Spectrum1D(flux=flux * u.Unit(start_unit), wcs=w, meta=wcs_dict)
cubeviz_helper.load_data(cube, data_label="test")

uc_plg = cubeviz_helper.plugins['Unit Conversion']
sv = cubeviz_helper.viewers['spectrum-viewer']
sv.set_limits(x_min=5e-7, x_max=6e-7, y_min=100, y_max=105)

uc_plg.flux_unit = end_unit
uc_plg.spectral_y_type = end_spectral_y_type

new_limits = sv._obj.get_limits()
assert np.allclose(new_limits, expected_limits)
Original file line number Diff line number Diff line change
Expand Up @@ -362,9 +362,6 @@ def _handle_spectral_y_unit(self, *args):
pass
else:
self.spectrum_viewer.set_plot_axes()
# until we can have upstream automatic limit updating on change
# in display units with equivalencies, we'll reset the limits
self.spectrum_viewer.reset_limits()

# broadcast that there has been a change in the spectrum viewer y axis,
self.hub.broadcast(GlobalDisplayUnitChanged('spectral_y',
Expand Down
34 changes: 33 additions & 1 deletion jdaviz/core/freezable_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,12 @@
from echo import delay_callback, CallbackProperty
import numpy as np

from astropy import units as u
from glue.viewers.profile.state import ProfileViewerState
from glue_jupyter.bqplot.image.state import BqplotImageViewerState
from glue.viewers.matplotlib.state import DeferredDrawCallbackProperty as DDCProperty

from jdaviz.utils import get_reference_image_data
from jdaviz.utils import get_reference_image_data, flux_conversion

__all__ = ['FreezableState', 'FreezableProfileViewerState', 'FreezableBqplotImageViewerState']

Expand Down Expand Up @@ -53,6 +54,37 @@ def _reset_x_limits(self, *event):
self.x_min = x_min
self.x_max = x_max

def _convert_units_y_limits(self, old_unit, new_unit):
# override glue's _convert_units_y_limits to account
# for equivalencies. This converts all four corners
# of the limits to set new limits that contain those
# same corners

if old_unit != new_unit:
if old_unit is None or new_unit is None:
self._reset_y_limits()
return
x_corners = np.array([self.x_min, self.x_min, self.x_max, self.x_max])
y_corners = np.array([self.y_min, self.y_max, self.y_min, self.y_max])
spectral_axis = x_corners * u.Unit(self.x_display_unit)

# NOTE: this uses the scale factor from the first found layer. We may want to
# generalize this to iterate over all scale factors if/when we support multiple
# flux cubes (with potentially different pixel scale factors).
for layer in self.layers:
if psc := getattr(layer.layer, 'meta', {}).get('_pixel_scale_factor', None): # noqa
spectral_axis.info.meta = {'_pixel_scale_factor',
psc}
break
else:
spectral_axis.info.meta = {}

y_corners_new = flux_conversion(y_corners, old_unit, new_unit, spectral_axis=spectral_axis) # noqa
Copy link
Contributor

Choose a reason for hiding this comment

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

will this handle equivalencies needed for surface brightness > surface brightness conversions, that may need pixel scale factor?

Copy link
Member Author

Choose a reason for hiding this comment

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

yes, it should

Copy link
Contributor

Choose a reason for hiding this comment

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

where does it get that from if there is no spectrum.meta or additional equivalency?

Copy link
Member Author

Choose a reason for hiding this comment

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

the for-loop directly above this (it does rely on layer.layer.meta, but that should be set by spectral extraction)

Copy link
Contributor

@gibsongreen gibsongreen Dec 13, 2024

Choose a reason for hiding this comment

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

On L74 there is a loop to find the scale factor if it exists, and if so, to attach it to the metadata of the spectral_axis. Then there is a small modification to flux_conversion to account for passing the scale factor to spectral_axis instead of spec(for the case that spec nor eqv is not passed).

Copy link
Contributor

Choose a reason for hiding this comment

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

thanks, i see what's happening now!


with delay_callback(self, 'y_min', 'y_max'):
self.y_min = np.nanmin(y_corners_new)
self.y_max = np.nanmax(y_corners_new)


class FreezableBqplotImageViewerState(BqplotImageViewerState, FreezableState):
linked_by_wcs = False
Expand Down
10 changes: 5 additions & 5 deletions jdaviz/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ def test_spec_sb_flux_conversion():
targ_units=(u.ph / (u.s * u.cm**2 * u.Hz * solid_angle)), # noqa
eqv=eqv,
spec_unit=spec_unit,
image_data=None
indirect_needs_spec_axis=None
)
assert_allclose(returned_values, target_values)
assert (return_units == expected_units)
Expand All @@ -71,7 +71,7 @@ def test_spec_sb_flux_conversion():
values=values,
orig_units=(u.ph / (u.Angstrom * u.s * u.cm**2 * u.sr)), # noqa
targ_units=(u.MJy), eqv=eqv,
spec_unit=spec_unit, image_data=None
spec_unit=spec_unit, indirect_needs_spec_axis=None # noqa
)
assert_allclose(returned_values, target_values)
assert (return_units == expected_units)
Expand All @@ -83,7 +83,7 @@ def test_spec_sb_flux_conversion():
returned_values, return_units, unit_flag = _indirect_conversion(
values=values, orig_units=(u.Jy/u.sr),
targ_units=(u.MJy), eqv=eqv,
spec_unit=spec_unit, image_data=None
spec_unit=spec_unit, indirect_needs_spec_axis=None # noqa
)
assert_allclose(returned_values, target_values)
assert (return_units == expected_units)
Expand All @@ -95,7 +95,7 @@ def test_spec_sb_flux_conversion():
returned_values, return_units = _indirect_conversion(
values=1, orig_units=(u.MJy/u.sr),
targ_units=(u.erg / (u.s * u.cm**2 * u.Hz * u.sr)),
eqv=eqv, spec_unit=None, image_data=True
eqv=eqv, spec_unit=None, indirect_needs_spec_axis=True
)
assert_allclose(returned_values, target_value)
assert return_units == expected_units
Expand All @@ -105,7 +105,7 @@ def test_spec_sb_flux_conversion():
expected_units = (u.MJy / u.sr)
returned_values, return_units = _indirect_conversion(
values=10, orig_units=(u.MJy/u.sr), targ_units=(u.Jy/u.sr),
eqv=eqv, spec_unit=None, image_data=True
eqv=eqv, spec_unit=None, indirect_needs_spec_axis=True
)
assert_allclose(returned_values, target_value)
assert return_units == expected_units
Expand Down
44 changes: 30 additions & 14 deletions jdaviz/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,7 @@ def indirect_units():
return units


def flux_conversion(values, original_units, target_units, spec=None, eqv=None, slice=None):
def flux_conversion(values, original_units, target_units, spec=None, eqv=None, spectral_axis=None):
"""
Convert flux or surface brightness values from original units to target units.

Expand All @@ -367,8 +367,8 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
eqv : list of :ref:`astropy:unit_equivalencies`, optional
A list of Astropy equivalencies necessary for complex unit conversions/translations.

slice : `astropy.units.Quantity`, optional
The current slice of a data cube, with units. Necessary for complex unit
spectral_axis : `astropy.units.Quantity`, optional
Provide spectral axis associated with values array. Necessary for simple and complex unit
conversions/translations that require spectral density equivalencies.

Returns
Expand All @@ -383,7 +383,7 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
# 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
image_data = False
indirect_needs_spec_axis = False
if spec:
if not np.isscalar(values) and len(values) == 2:
spectral_values = spec.spectral_axis[0]
Expand All @@ -398,10 +398,18 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
eqv = u.spectral_density(spectral_values)
else:
eqv = u.spectral_density(spec.spectral_axis[0])
elif slice is not None and eqv:
image_data = True
# Need this to convert Flux to Flux for complex conversions/translations of cube image data
eqv += u.spectral_density(slice)
elif spectral_axis is not None:
if not isinstance(spectral_axis, u.Quantity):
raise TypeError("spectral_axis must be an instance of astropy.units.Quantity.")

if eqv:
indirect_needs_spec_axis = True
# Needed to convert Flux to Flux for complex conversion/translation of non-spectrum data
eqv += u.spectral_density(spectral_axis)
elif not np.isscalar(values) and len(values) != spectral_axis.size:
raise ValueError(f"Size mismatch: values size : ({values.size}) must match spectral_axis size : ({spectral_axis.size}).") # noqa
else:
eqv = u.spectral_density(spectral_axis)

orig_units = u.Unit(original_units)
targ_units = u.Unit(target_units)
Expand All @@ -410,7 +418,9 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
solid_angle_in_targ = check_if_unit_is_per_solid_angle(targ_units, return_unit=True)

# Ensure a spectrum passed through Spectral Extraction plugin
if (((spec and ('_pixel_scale_factor' in spec.meta))) and
if ((((spec and ('_pixel_scale_factor' in spec.meta))) or
(spectral_axis is not None and ('_pixel_scale_factor' in spectral_axis.info.meta)))
and
(((solid_angle_in_orig) and (not solid_angle_in_targ)) or
((not solid_angle_in_orig) and (solid_angle_in_targ)))):
# Data item in data collection does not update from conversion/translation.
Expand All @@ -419,7 +429,13 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
n_values = len(values)

# Make sure they are float (can be Quantity).
fac = spec.meta['_pixel_scale_factor']
if spec:
fac = spec.meta['_pixel_scale_factor']
elif spectral_axis is not None:
# Quantity.info.meta gets casted to set
fac = next((item for item in spectral_axis.info.meta if isinstance(item, float)), None)
spec_unit = orig_units

if isinstance(fac, Quantity):
fac = fac.value

Expand Down Expand Up @@ -454,10 +470,10 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s
elif selected_unit_updated == 'orig':
orig_units = updated_units

elif image_data:
elif indirect_needs_spec_axis:
values, orig_units = _indirect_conversion(
values=values, orig_units=orig_units, targ_units=targ_units,
eqv=eqv, image_data=image_data
eqv=eqv, indirect_needs_spec_axis=indirect_needs_spec_axis
)
elif solid_angle_in_orig == solid_angle_in_targ == PIX2:
# in the case where we have 2 SBs per solid pixel that need
Expand All @@ -473,7 +489,7 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s


def _indirect_conversion(values, orig_units, targ_units, eqv,
spec_unit=None, image_data=None):
spec_unit=None, indirect_needs_spec_axis=None):

# Note: is there a way we could write this to not require 'spec_unit'? It
# seems like it falls back on this to get a solid angle unit, but can we
Expand Down Expand Up @@ -502,7 +518,7 @@ def _indirect_conversion(values, orig_units, targ_units, eqv,

return values, targ_units, 'targ'

elif image_data or (spec_unit and solid_angle_in_spec):
elif indirect_needs_spec_axis or (spec_unit and solid_angle_in_spec):
if not solid_angle_in_targ:
targ_units /= solid_angle_in_spec
if ((u.Unit(targ_units) in indirect_units()) or
Expand Down
Loading