diff --git a/CHANGES.rst b/CHANGES.rst index 2d6fb7cc6d..145c8f02c9 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -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] + #3139, #3149, #3155, #3178, #3185, #3187, #3190, #3156] - Plugin tray is now open by default. [#2892] diff --git a/jdaviz/app.py b/jdaviz/app.py index 6df7f76217..17c0981c2d 100644 --- a/jdaviz/app.py +++ b/jdaviz/app.py @@ -78,15 +78,27 @@ def equivalent_units(self, data, cid, units): 'Jy', 'mJy', 'uJy', 'MJy', 'W / (Hz m2)', 'eV / (Hz s m2)', 'erg / (Hz s cm2)', 'erg / (Angstrom s cm2)', - 'ph / (Angstrom s cm2)', 'ph / (Hz s cm2)' + 'ph / (Angstrom s cm2)', 'ph / (Hz s cm2)', + 'ct' ] + [ 'Jy / sr', 'mJy / sr', 'uJy / sr', 'MJy / sr', 'W / (Hz sr m2)', 'eV / (Hz s sr m2)', 'erg / (s sr cm2)', 'erg / (Hz s sr cm2)', 'erg / (Angstrom s sr cm2)', - 'ph / (Angstrom s sr cm2)', 'ph / (Hz s sr cm2)' - ]) + 'ph / (Angstrom s sr cm2)', 'ph / (Hz s sr cm2)', + 'ct / sr' + ] + + [ + 'Jy / pix2', 'mJy / pix2', 'uJy / pix2', 'MJy / pix2', + 'W / (Hz m2 pix2)', 'eV / (Hz s m2 pix2)', + 'erg / (s cm2 pix2)', 'erg / (Hz s cm2 pix2)', + 'erg / (Angstrom s cm2 pix2)', + 'ph / (Angstrom s cm2 pix2)', 'ph / (Hz s cm2 pix2)', + 'ct / pix2' + ] + + ) else: # spectral axis # prefer Hz over Bq and um over micron exclude = {'Bq', 'micron'} @@ -1294,16 +1306,38 @@ def _get_display_unit(self, axis): sv_y_unit = sv.data()[0].flux.unit if axis == 'spectral_y': return sv_y_unit - elif axis == 'flux': - if check_if_unit_is_per_solid_angle(sv_y_unit): - # TODO: this will need updating once solid-angle unit can be non-steradian - return sv_y_unit * u.sr + # since this is where we're falling back on native units (UC plugin might not exist) + # first check the spectrum viewer y axis for any solid angle unit (i think that it + # will ALWAYS be in flux, but just to be sure). If no solid angle unit is found, + # check the flux viewer for surface brightness units + sv_y_angle_unit = check_if_unit_is_per_solid_angle(sv_y_unit, return_unit=True) + + # check flux viewer if none in spectral viewer + fv_angle_unit = None + if not sv_y_angle_unit: + if hasattr(self._jdaviz_helper, '_default_flux_viewer_reference_name'): + vname = self._jdaviz_helper._default_flux_viewer_reference_name + fv = self.get_viewer(vname) + fv_unit = fv.data()[0].get_object().flux.unit + fv_angle_unit = check_if_unit_is_per_solid_angle(fv_unit, + return_unit=True) + else: + # mosviz, not sure what to do here but can't access flux + # viewer the same way. once we force the UC plugin to + # exist this won't matter anyway because units can be + # acessed from the plugin directly. assume u.sr for now + fv_angle_unit = u.sr + + solid_angle_unit = sv_y_angle_unit or fv_angle_unit + + if axis == 'flux': + if sv_y_angle_unit: + return sv_y_unit * solid_angle_unit return sv_y_unit - else: - # surface_brightness - if check_if_unit_is_per_solid_angle(sv_y_unit): + elif axis == 'sb': + if sv_y_angle_unit: return sv_y_unit - return sv_y_unit / u.sr + return sv_y_unit / solid_angle_unit elif axis == 'temporal': # No unit for ramp's time (group/resultant) axis: return None diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py index 9cdd826847..7748a3511d 100644 --- a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py +++ b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py @@ -20,7 +20,6 @@ SpectralContinuumMixin, skip_if_no_updates_since_last_active, with_spinner) -from jdaviz.core.validunits import check_if_unit_is_per_solid_angle from jdaviz.core.user_api import PluginUserApi __all__ = ['MomentMap'] @@ -111,8 +110,10 @@ def __init__(self, *args, **kwargs): self.output_unit = SelectPluginComponent(self, items='output_unit_items', selected='output_unit_selected', - manual_options=['Flux', 'Spectral Unit', - 'Velocity', 'Velocity^N']) + manual_options=['Surface Brightness', + 'Spectral Unit', + 'Velocity', + 'Velocity^N']) self.dataset.add_filter('is_cube') self.add_results.viewer.filters = ['is_image_viewer'] @@ -174,16 +175,11 @@ def _set_data_units(self, event={}): self.output_unit_selected = moment_unit_options[unit_options_index][0] self.send_state("output_unit_selected") - # either 'Flux' or 'Surface Brightness' - orig_spectral_y_type = self.output_unit_items[0]['label'] - - unit_dict = {orig_spectral_y_type: "", + unit_dict = {"Surface Brightness": "", "Spectral Unit": "", "Velocity": "km/s", "Velocity^N": f"km{self.n_moment}/s{self.n_moment}"} - sb_or_flux_label = None - if self.dataset_selected != "": # Spectral axis is first in this list data = self.app.data_collection[self.dataset_selected] @@ -196,37 +192,11 @@ def _set_data_units(self, event={}): sunit = "" self.dataset_spectral_unit = sunit unit_dict["Spectral Unit"] = sunit - - # get flux/SB units - if self.spectrum_viewer and hasattr(self.spectrum_viewer.state, 'y_display_unit'): - if self.spectrum_viewer.state.y_display_unit is not None: - unit_dict[orig_spectral_y_type] = self.spectrum_viewer.state.y_display_unit - else: - # spectrum_viewer.state will only have x/y_display_unit if unit conversion has - # been done if not, get default flux units which should be the units displayed - unit_dict[orig_spectral_y_type] = data.get_component('flux').units - else: - # spectrum_viewer.state will only have x/y_display_unit if unit conversion has - # been done if not, get default flux units which should be the units displayed - unit_dict[orig_spectral_y_type] = data.get_component('flux').units - - # figure out if label should say 'Flux' or 'Surface Brightness' - sb_or_flux_label = "Flux" - is_unit_solid_angle = check_if_unit_is_per_solid_angle(unit_dict[orig_spectral_y_type]) - if is_unit_solid_angle is True: - sb_or_flux_label = "Surface Brightness" + unit_dict["Surface Brightness"] = self.app._get_display_unit('sb') # Update units in selection item dictionary for item in self.output_unit_items: item["unit_str"] = unit_dict[item["label"]] - if item["label"] in ["Flux", "Surface Brightness"] and sb_or_flux_label: - # change unit label to reflect if unit is flux or SB - item["label"] = sb_or_flux_label - - if self.dataset_selected != "": - # output_unit_selected might not match (potentially) changed flux/SB label - if self.output_unit_selected in ['Flux', 'Surface Brightness']: - self.output_unit_selected = sb_or_flux_label # Filter what we want based on n_moment if self.n_moment == 0: @@ -343,18 +313,18 @@ def calculate_moment(self, add_data=True): # convert units for moment 0, which is the only currently supported # moment for using converted units. if n_moment == 0: - # get flux/SB units - if self.spectrum_viewer and hasattr(self.spectrum_viewer.state, 'y_display_unit'): - if self.spectrum_viewer.state.y_display_unit is not None: - flux_sb_unit = self.spectrum_viewer.state.y_display_unit - else: - flux_sb_unit = data.get_component('flux').units - else: - flux_sb_unit = data.get_component('flux').units + # get display units for moment 0 based on unit conversion plugin selection + # the 0th item of this dictionary is the surface brightness unit, kept + # up to date with choices from the UC plugin + moment_0_display_unit = self.output_unit_items[0]['unit_str'] + + # convert unit string to Unit so moment map data can be converted + # `display_unit` is the data unit (surface brighntess) translated + # to the choice of selected flux and angle unit from the UC plugin. + display_unit = u.Unit(moment_0_display_unit) + + moment_new_unit = display_unit * self.app._get_display_unit('spectral') # noqa: E501 - # convert unit string to u.Unit so moment map data can be converted - spectral_y_display_unit = u.Unit(flux_sb_unit) - moment_new_unit = spectral_y_display_unit * self.spectrum_viewer.state.x_display_unit # noqa: E501 self.moment = self.moment.to(moment_new_unit) # Reattach the WCS so we can load the result diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py index 90b35bb3db..fa6941b81a 100644 --- a/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py +++ b/jdaviz/configs/cubeviz/plugins/moment_maps/tests/test_moment_maps.py @@ -12,10 +12,21 @@ from numpy.testing import assert_allclose -def test_user_api(cubeviz_helper, spectrum1d_cube): +@pytest.mark.parametrize("cube_type", ["Surface Brightness", "Flux"]) +def test_user_api(cubeviz_helper, spectrum1d_cube, spectrum1d_cube_sb_unit, cube_type): + + # test is parameterize to test a cube that is in Jy / sr (Surface Brightness) + # as well as Jy (Flux), to test that flux cubes, which are converted in the + # parser to flux / pix^2 surface brightness cubes, both work correctly. + + if cube_type == 'Surface Brightness': + cube = spectrum1d_cube_sb_unit + elif cube_type == 'Flux': + cube = spectrum1d_cube + with warnings.catch_warnings(): warnings.filterwarnings("ignore", message="No observer defined on WCS.*") - cubeviz_helper.load_data(spectrum1d_cube, data_label='test') + cubeviz_helper.load_data(cube, data_label='test') mm = cubeviz_helper.plugins['Moment Maps'] assert not mm._obj.continuum_marks['center'].visible @@ -43,20 +54,33 @@ def test_user_api(cubeviz_helper, spectrum1d_cube): mm.n_moment = 1 with pytest.raises(ValueError, match="not one of"): mm._obj.output_unit_selected = "Bad Unit" - mm._obj.output_unit_selected = "Flux" + mm._obj.output_unit_selected = "Surface Brightness" with pytest.raises(ValueError, match="units must be in"): mm.calculate_moment() -def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path): +@pytest.mark.parametrize("cube_type", ["Surface Brightness", "Flux"]) +def test_moment_calculation(cubeviz_helper, spectrum1d_cube, + spectrum1d_cube_sb_unit, cube_type, tmp_path): moment_unit = "Jy m" moment_value_str = "+6.40166e-10" + if cube_type == 'Surface Brightness': + moment_unit += " / sr" + cube = spectrum1d_cube_sb_unit + cube_unit = cube.unit.to_string() + + elif cube_type == 'Flux': + moment_unit += " / pix2" + cube = spectrum1d_cube + cube_unit = cube.unit.to_string() + " / pix2" # cube in Jy will become cube in Jy / pix2 + dc = cubeviz_helper.app.data_collection with warnings.catch_warnings(): warnings.filterwarnings("ignore", message="No observer defined on WCS.*") - cubeviz_helper.load_data(spectrum1d_cube, data_label='test') + cubeviz_helper.load_data(cube, data_label='test') + flux_viewer = cubeviz_helper.app.get_viewer(cubeviz_helper._default_flux_viewer_reference_name) # Since we are not really displaying, need this to trigger GUI stuff. @@ -68,7 +92,9 @@ def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path): mm.n_moment = 0 # Collapsed sum, will get back 2D spatial image assert mm._obj.results_label == 'moment 0' - assert mm.output_unit == "Flux" + + # result is always a SB unit - if flux cube loaded, per pix2 + assert mm.output_unit == 'Surface Brightness' mm._obj.add_results.viewer.selected = cubeviz_helper._default_uncert_viewer_reference_name mm._obj.vue_calculate_moment() @@ -88,13 +114,13 @@ def test_moment_calculation(cubeviz_helper, spectrum1d_cube, tmp_path): assert mm._obj.results_label == 'moment 0' assert mm._obj.results_label_overwrite is True - # Make sure coordinate display works + # Make sure coordinate display works in flux viewer (loaded data, not the moment map) label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 0, 'y': 0}}) assert flux_viewer.state.slices == (0, 0, 1) # Slice 0 has 8 pixels, this is Slice 1 - assert label_mouseover.as_text() == ("Pixel x=00.0 y=00.0 Value +8.00000e+00 Jy", + assert label_mouseover.as_text() == (f"Pixel x=00.0 y=00.0 Value +8.00000e+00 {cube_unit}", "World 13h39m59.9731s +27d00m00.3600s (ICRS)", "204.9998877673 27.0001000000 (deg)") @@ -290,7 +316,7 @@ def test_momentmap_nirspec_prism(cubeviz_helper, tmp_path): (sky_cube.ra.deg, sky_cube.dec.deg)) -def test_correct_output_flux_or_sb_units(cubeviz_helper, spectrum1d_cube_custom_fluxunit): +def test_correct_output_spectral_y_units(cubeviz_helper, spectrum1d_cube_custom_fluxunit): moment_unit = "Jy m / sr" diff --git a/jdaviz/configs/cubeviz/plugins/parsers.py b/jdaviz/configs/cubeviz/plugins/parsers.py index d9b2484a11..95410f239b 100644 --- a/jdaviz/configs/cubeviz/plugins/parsers.py +++ b/jdaviz/configs/cubeviz/plugins/parsers.py @@ -8,9 +8,10 @@ from astropy.nddata import StdDevUncertainty from astropy.time import Time from astropy.wcs import WCS -from specutils import Spectrum1D +from specutils import Spectrum1D, SpectralAxis from jdaviz.core.registries import data_parser_registry +from jdaviz.core.validunits import check_if_unit_is_per_solid_angle from jdaviz.utils import standardize_metadata, PRIHDR_KEY, download_uri_to_path @@ -177,15 +178,27 @@ def _get_celestial_wcs(wcs): return wcs.celestial if hasattr(wcs, 'celestial') else None -def _return_spectrum_with_correct_units(flux, wcs, metadata, data_type, target_wave_unit=None, - hdulist=None, uncertainty=None, mask=None): +def _return_spectrum_with_correct_units(flux, wcs, metadata, data_type, + target_wave_unit=None, hdulist=None, + uncertainty=None, mask=None, apply_pix2=False): """Upstream issue of WCS not using the correct units for data must be fixed here. - Issue: https://github.com/astropy/astropy/issues/3658 + Issue: https://github.com/astropy/astropy/issues/3658. + + Also converts flux units to flux/pix2 solid angle units, if `flux` is not a surface + brightness and `apply_pix2` is True. """ with warnings.catch_warnings(): warnings.filterwarnings( 'ignore', message='Input WCS indicates that the spectral axis is not last', category=UserWarning) + + # convert flux and uncertainty to per-pix2 if input is not a surface brightness + if apply_pix2: + if not check_if_unit_is_per_solid_angle(flux.unit): + flux = flux / (u.pix * u.pix) + if uncertainty is not None: + uncertainty = uncertainty / (u.pix * u.pix) + sc = Spectrum1D(flux=flux, wcs=wcs, uncertainty=uncertainty, mask=mask) if target_wave_unit is None and hdulist is not None: @@ -281,7 +294,9 @@ def _parse_hdulist(app, hdulist, file_name=None, # to sky regions, where the parent data of the subset might have dropped spatial WCS info metadata['_orig_spatial_wcs'] = _get_celestial_wcs(wcs) - sc = _return_spectrum_with_correct_units(flux, wcs, metadata, data_type, hdulist=hdulist) + apply_pix2 = data_type in ['flux', 'uncert'] + sc = _return_spectrum_with_correct_units(flux, wcs, metadata, data_type, + hdulist=hdulist, apply_pix2=apply_pix2) app.add_data(sc, data_label) @@ -295,7 +310,8 @@ def _parse_hdulist(app, hdulist, file_name=None, else: # flux # Forced wave unit conversion made it lose stuff, so re-add - app.data_collection[data_label].get_component("flux").units = flux_unit + # also re-get unit from sc in case a factor of pix2 was applied + app.data_collection[data_label].get_component("flux").units = sc.unit # Add flux to top left image viewer app.add_data_to_viewer(flux_viewer_reference_name, data_label) app._jdaviz_helper._loaded_flux_cube = app.data_collection[data_label] @@ -341,7 +357,6 @@ def _parse_jwst_s3d(app, hdulist, data_label, ext='SCI', app.add_data(data, data_label, parent=parent) # get glue data and update if DQ: - if ext == 'DQ': # prevent circular import: from jdaviz.configs.imviz.plugins.parsers import prep_data_layer_as_dq @@ -448,6 +463,11 @@ def _parse_spectrum1d_3d(app, file_obj, data_label=None, s1d = Spectrum1D(flux=flux, wcs=file_obj.wcs, meta=meta) + # convert data loaded in flux units to a per-square-pixel surface + # brightness unit (e.g Jy to Jy/pix**2) + if (attr != "mask") and (not check_if_unit_is_per_solid_angle(flux.unit)): + s1d = convert_spectrum1d_from_flux_to_flux_per_pixel(s1d) + cur_data_label = app.return_data_label(data_label, attr.upper()) app.add_data(s1d, cur_data_label) @@ -474,6 +494,11 @@ def _parse_spectrum1d(app, file_obj, data_label=None, spectrum_viewer_reference_ # TODO: glue-astronomy translators only look at the flux property of # specutils Spectrum1D objects. Fix to support uncertainties and masks. + # convert data loaded in flux units to a per-square-pixel surface + # brightness unit (e.g Jy to Jy/pix**2) + if not check_if_unit_is_per_solid_angle(file_obj.flux.unit): + file_obj = convert_spectrum1d_from_flux_to_flux_per_pixel(file_obj) + app.add_data(file_obj, data_label) app.add_data_to_viewer(spectrum_viewer_reference_name, data_label) @@ -496,6 +521,12 @@ def _parse_ndarray(app, file_obj, data_label=None, data_type=None, meta = standardize_metadata({'_orig_spatial_wcs': None}) s3d = Spectrum1D(flux=flux, meta=meta) + + # convert data loaded in flux units to a per-square-pixel surface + # brightness unit (e.g Jy to Jy/pix**2) + if not check_if_unit_is_per_solid_angle(s3d.unit): + file_obj = convert_spectrum1d_from_flux_to_flux_per_pixel(s3d) + app.add_data(s3d, data_label) if data_type == 'flux': @@ -522,6 +553,11 @@ def _parse_gif(app, file_obj, data_label=None, flux_viewer_reference_name=None): meta = {'filename': file_name, '_orig_spatial_wcs': None} s3d = Spectrum1D(flux=flux * u.count, meta=standardize_metadata(meta)) + # convert data loaded in flux units to a per-square-pixel surface + # brightness unit (e.g Jy to Jy/pix**2) + if not check_if_unit_is_per_solid_angle(s3d): + file_obj = convert_spectrum1d_from_flux_to_flux_per_pixel(s3d) + app.add_data(s3d, data_label) app.add_data_to_viewer(flux_viewer_reference_name, data_label) @@ -539,3 +575,66 @@ def _get_data_type_by_hdu(hdu): else: data_type = '' return data_type + + +def convert_spectrum1d_from_flux_to_flux_per_pixel(spectrum): + """ + Converts a Spectrum1D object's flux units to flux per square pixel. + + This function takes a `specutils.Spectrum1D` object with flux units and converts the + flux (and optionally, uncertainty) to a surface brightness per square pixel + (e.g., from Jy to Jy/pix**2). This is done by updating the units of spectrum.flux + and (if present) spectrum.uncertainty, and creating a new `specutils.Spectrum1D` + object with the modified flux and uncertainty. + + Parameters + ---------- + spectrum : Spectrum1D + A `specutils.Spectrum1D` object containing flux data, which is assumed to be in + flux units without any angular component in the denominator. + + Returns + ------- + Spectrum1D + A new `specutils.Spectrum1D` object with flux and uncertainty (if present) + converted to units of flux per square pixel. + """ + + # convert flux, which is always populated + flux = getattr(spectrum, 'flux') + flux = flux / (u.pix * u.pix) + + # and uncerts, if present + uncerts = getattr(spectrum, 'uncertainty') + if uncerts is not None: + old_uncerts = uncerts.represent_as(StdDevUncertainty) # enforce common uncert type. + uncerts = old_uncerts.quantity / (u.pix * u.pix) + + # create a new spectrum 1d with all the info from the input spectrum 1d, + # and the flux / uncerts converted from flux to SB per square pixel + + # if there is a spectral axis that is a SpectralAxis, you cant also set + # redshift or radial_velocity + spectral_axis = getattr(spectrum, 'spectral_axis', None) + if spectral_axis is not None: + if isinstance(spectral_axis, SpectralAxis): + redshift = None + radial_velocity = None + else: + redshift = spectrum.redshift + radial_velocity = spectrum.radial_velocity + + # initialize new spectrum1d with new flux, uncerts, and all other init parameters + # from old input spectrum as well as any 'meta'. any more missing information + # not in init signature that might be present in `spectrum`? + new_spec1d = Spectrum1D(flux=flux, uncertainty=uncerts, + spectral_axis=spectrum.spectral_axis, + mask=spectrum.mask, + wcs=spectrum.wcs, + velocity_convention=spectrum.velocity_convention, + rest_value=spectrum.rest_value, redshift=redshift, + radial_velocity=radial_velocity, + bin_specification=getattr(spectrum, 'bin_specification', None), + meta=spectrum.meta) + + return new_spec1d diff --git a/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py b/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py index b649e9f06a..a01101dfdb 100644 --- a/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py +++ b/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py @@ -3,7 +3,6 @@ import numpy as np import astropy -import astropy.units as u from astropy.nddata import ( NDDataArray, StdDevUncertainty ) @@ -24,6 +23,7 @@ skip_if_no_updates_since_last_active, with_spinner, with_temp_disable) from jdaviz.core.user_api import PluginUserApi +from jdaviz.core.validunits import check_if_unit_is_per_solid_angle from jdaviz.configs.cubeviz.plugins.parsers import _return_spectrum_with_correct_units from jdaviz.configs.cubeviz.plugins.viewers import WithSliceIndicator @@ -506,9 +506,16 @@ def _extract_from_aperture(self, cube, uncert_cube, aperture, collapsed_nddata = getattr(nddata_reshaped, selected_func)( axis=self.spatial_axes, **kwargs ) # returns an NDDataArray - # Remove per steradian denominator - if astropy.units.sr in collapsed_nddata.unit.bases: - aperture_area = self.cube.meta.get('PIXAR_SR', 1.0) * u.sr + + # Remove per solid angle denominator to turn sb into flux + sq_angle_unit = check_if_unit_is_per_solid_angle(collapsed_nddata.unit, + return_unit=True) + if sq_angle_unit is not None: + # convert aperture area in steradians to the selected square angle unit + # NOTE: just forcing these units for now!! this is in steradians and + # needs to be converted to the selected square angle unit but for now just + # force to correct units + aperture_area = self.cube.meta.get('PIXAR_SR', 1.0) * sq_angle_unit collapsed_nddata = collapsed_nddata.multiply(aperture_area, propagate_uncertainties=True) else: diff --git a/jdaviz/configs/cubeviz/plugins/spectral_extraction/tests/test_spectral_extraction.py b/jdaviz/configs/cubeviz/plugins/spectral_extraction/tests/test_spectral_extraction.py index 0cfc4044d6..084bbca5bf 100644 --- a/jdaviz/configs/cubeviz/plugins/spectral_extraction/tests/test_spectral_extraction.py +++ b/jdaviz/configs/cubeviz/plugins/spectral_extraction/tests/test_spectral_extraction.py @@ -35,6 +35,11 @@ def test_version_after_nddata_update(cubeviz_helper, spectrum1d_cube_with_uncert # Axes 0, 1 are the spatial ones. collapsed_cube_nddata = spectral_cube.sum(axis=(0, 1)) # return NDDataArray + # when loaded into app, cubes loaded in flux are converted to per-pixel-squared + # surface brightness, so multiply by pix**2 to compare to NDData, if input + # cube was in flux + collapsed_cube_nddata = collapsed_cube_nddata * (u.pix ** 2) + # Collapse the spectral cube using the methods in jdaviz: collapsed_cube_s1d = plg.extract(add_data=False) # returns Spectrum1D diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_cubeviz_aperphot.py b/jdaviz/configs/cubeviz/plugins/tests/test_cubeviz_aperphot.py index 8704cd0c62..43b170c3e6 100644 --- a/jdaviz/configs/cubeviz/plugins/tests/test_cubeviz_aperphot.py +++ b/jdaviz/configs/cubeviz/plugins/tests/test_cubeviz_aperphot.py @@ -10,7 +10,8 @@ def test_cubeviz_aperphot_cube_orig_flux(cubeviz_helper, image_cube_hdu_obj_microns): cubeviz_helper.load_data(image_cube_hdu_obj_microns, data_label="test") - flux_unit = u.Unit("1E-17 erg*s^-1*cm^-2*Angstrom^-1") + flux_unit = u.Unit("1E-17 erg*s^-1*cm^-2*Angstrom^-1*pix^-2") # actually a sb + solid_angle_unit = u.pix * u.pix aper = RectanglePixelRegion(center=PixCoord(x=1, y=2), width=3, height=5) cubeviz_helper.load_regions(aper) @@ -33,8 +34,11 @@ def test_cubeviz_aperphot_cube_orig_flux(cubeviz_helper, image_cube_hdu_obj_micr sky = row["sky_center"] assert_allclose(sky.ra.deg, 205.43985906934287) assert_allclose(sky.dec.deg, 27.003490103642033) - assert_allclose(row["sum"], 75 * flux_unit) # 3 (w) x 5 (h) x 5 (v) - assert_allclose(row["sum_aper_area"], 15 * (u.pix * u.pix)) # 3 (w) x 5 (h) + + # sum should be in flux ( have factor of pix^2 multiplied out of input unit) + assert_allclose(row["sum"], 75 * flux_unit * solid_angle_unit) # 3 (w) x 5 (h) x 5 (v) + + assert_allclose(row["sum_aper_area"], 15 * solid_angle_unit) # 3 (w) x 5 (h) assert_allclose(row["mean"], 5 * flux_unit) assert_quantity_allclose(row["slice_wave"], 4.894499866699333 * u.um) @@ -50,8 +54,11 @@ def test_cubeviz_aperphot_cube_orig_flux(cubeviz_helper, image_cube_hdu_obj_micr sky = row["sky_center"] assert_allclose(sky.ra.deg, 205.43985906934287) assert_allclose(sky.dec.deg, 27.003490103642033) - assert_allclose(row["sum"], 15 * flux_unit) # 3 (w) x 5 (h) x 1 (v) - assert_allclose(row["sum_aper_area"], 15 * (u.pix * u.pix)) # 3 (w) x 5 (h) + + # sum should be in flux ( have factor of pix^2 multiplied out of input unit) + assert_allclose(row["sum"], 15 * flux_unit * solid_angle_unit) # 3 (w) x 5 (h) x 1 (v) + + assert_allclose(row["sum_aper_area"], 15 * solid_angle_unit) # 3 (w) x 5 (h) assert_allclose(row["mean"], 1 * flux_unit) assert_quantity_allclose(row["slice_wave"], 4.8904998665093435 * u.um) @@ -75,8 +82,11 @@ def test_cubeviz_aperphot_cube_orig_flux(cubeviz_helper, image_cube_hdu_obj_micr sky = row["sky_center"] assert_allclose(sky.ra.deg, 205.43985906934287) assert_allclose(sky.dec.deg, 27.003490103642033) - assert_allclose(row["sum"], 540 * flux_unit) # 3 (w) x 5 (h) x 36 (v) - assert_allclose(row["sum_aper_area"], 15 * (u.pix * u.pix)) # 3 (w) x 5 (h) + + # sum should be in flux ( have factor of pix^2 multiplied out of input unit) + assert_allclose(row["sum"], 540 * flux_unit * solid_angle_unit) # 3 (w) x 5 (h) x 36 (v) + + assert_allclose(row["sum_aper_area"], 15 * solid_angle_unit) # 3 (w) x 5 (h) assert_allclose(row["mean"], 36 * flux_unit) assert np.isnan(row["slice_wave"]) @@ -88,7 +98,8 @@ def test_cubeviz_aperphot_cube_orig_flux(cubeviz_helper, image_cube_hdu_obj_micr def test_cubeviz_aperphot_generated_3d_gaussian_smooth(cubeviz_helper, image_cube_hdu_obj_microns): cubeviz_helper.load_data(image_cube_hdu_obj_microns, data_label="test") - flux_unit = u.Unit("1E-17 erg*s^-1*cm^-2*Angstrom^-1") + flux_unit = u.Unit("1E-17 erg*s^-1*cm^-2*Angstrom^-1*pix^-2") # actually a sb + solid_angle_unit = u.pix * u.pix gauss_plg = cubeviz_helper.plugins["Gaussian Smooth"]._obj gauss_plg.mode_selected = "Spatial" @@ -113,13 +124,85 @@ def test_cubeviz_aperphot_generated_3d_gaussian_smooth(cubeviz_helper, image_cub sky = row["sky_center"] assert_allclose(sky.ra.deg, 205.43985906934287) assert_allclose(sky.dec.deg, 27.003490103642033) - assert_allclose(row["sum"], 48.54973 * flux_unit) # 3 (w) x 5 (h) x <5 (v) - assert_allclose(row["sum_aper_area"], 15 * (u.pix * u.pix)) # 3 (w) x 5 (h) + + # sum should be in flux ( have factor of pix^2 multiplied out of input unit) + assert_allclose(row["sum"], 48.54973 * flux_unit * solid_angle_unit) # 3 (w) x 5 (h) x <5 (v) + + assert_allclose(row["sum_aper_area"], 15 * solid_angle_unit) # 3 (w) x 5 (h) assert_allclose(row["mean"], 3.236648941040039 * flux_unit) assert_quantity_allclose(row["slice_wave"], 4.894499866699333 * u.um) +@pytest.mark.parametrize("cube_unit", [u.MJy / u.sr, u.MJy, u.MJy / (u.pix*u.pix)]) +def test_cubeviz_aperphot_cube_sr_and_pix2(cubeviz_helper, + spectrum1d_cube_custom_fluxunit, + cube_unit): + # tests aperture photometry outputs between different inputs of flux (which + # should be converted to a surface brighntess after loading), flux per sr + # and flux per square pixel. the pixel area for per steradian cubes is + # set so the output values between units will be the same + + cube = spectrum1d_cube_custom_fluxunit(fluxunit=cube_unit) + cubeviz_helper.load_data(cube, data_label="test") + + aper = RectanglePixelRegion(center=PixCoord(x=3, y=1), width=1, height=1) + bg = RectanglePixelRegion(center=PixCoord(x=2, y=0), width=1, height=1) + cubeviz_helper.load_regions([aper, bg]) + + plg = cubeviz_helper.plugins["Aperture Photometry"]._obj + plg.dataset_selected = "test[FLUX]" + plg.aperture_selected = "Subset 1" + plg.background_selected = "Subset 2" + + # Check that the default flux scaling is present for MJy / sr cubes + if cube_unit == (u.MJy / u.sr): + assert_allclose(plg.flux_scaling, 0.003631) + else: + assert_allclose(plg.flux_scaling, 0.0) + + # if cube is MJy / sr, set pixel area to 1 sr / pix2 so + # we can directly compare outputs between per sr and per pixel cubes, which + # will give the same results with this scaling + if cube_unit == (u.MJy / u.sr): + assert_allclose(plg.pixel_area, 0.01) # check default + plg.pixel_area = 1 * (u.sr).to(u.arcsec*u.arcsec) + solid_angle_unit = u.sr + + else: + # for per pixel cubes, set flux scaling to default for MJy / sr cubes + # so we can directly compare. this shouldn't be populated automatically, + # which is checked above + plg.flux_scaling = 0.003631 + solid_angle_unit = u.pix * u.pix + cube_unit = u.MJy / solid_angle_unit # cube unit in app is now per pix2 + + plg.vue_do_aper_phot() + row = cubeviz_helper.get_aperture_photometry_results()[0] + + # Basically, we should recover the input rectangle here, minus background. + assert_allclose(row["xcenter"], 3 * u.pix) + assert_allclose(row["ycenter"], 1 * u.pix) + # (15 - 10) MJy/sr x 1 sr, will always be MJy since solid angle is multiplied out + assert_allclose(row["sum"], 5.0 * u.MJy) + + assert_allclose(row["sum_aper_area"], 1 * (u.pix * u.pix)) + + # we forced area to be one sr so MJy / sr and MJy / pix2 gave the same result + assert_allclose(row["pixarea_tot"], 1.0 * solid_angle_unit) + + # also forced flux_scaling to be the same for MJy / sr cubes, which get a + # default value populated, and MJy / pix2 which have a default of 0.0 + assert_allclose(row["aperture_sum_mag"], -7.847359 * u.mag) + + assert_allclose(row["mean"], 5 * (cube_unit)) + # TODO: check if slice plugin has value_unit set correctly + assert_quantity_allclose(row["slice_wave"], 0.46236 * u.um) + + def test_cubeviz_aperphot_cube_orig_flux_mjysr(cubeviz_helper, spectrum1d_cube_custom_fluxunit): + # this test is essentially the same as test_cubeviz_aperphot_cube_sr_and_pix2 but for a single + # surface brightness unit and without changing the pixel area to make outputs the same. it + # was requested in review to keep both tests cube = spectrum1d_cube_custom_fluxunit(fluxunit=u.MJy / u.sr) cubeviz_helper.load_data(cube, data_label="test") diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py b/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py index 8946cc283e..40aa0c8efd 100644 --- a/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py +++ b/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py @@ -45,7 +45,7 @@ def test_fits_image_hdu_with_microns(image_cube_hdu_obj_microns, cubeviz_helper) label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 0, 'y': 0}}) - flux_unit_str = "erg / (Angstrom s cm2)" + flux_unit_str = "erg / (Angstrom s cm2 pix2)" assert label_mouseover.as_text() == (f'Pixel x=00.0 y=00.0 Value +5.00000e+00 1e-17 {flux_unit_str}', # noqa 'World 13h41m45.5759s +27d00m12.3044s (ICRS)', '205.4398995981 27.0034178810 (deg)') # noqa @@ -98,7 +98,7 @@ def test_fits_image_hdu_parse_from_file(tmpdir, image_cube_hdu_obj, cubeviz_help label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 0, 'y': 0}}) - flux_unit_str = "erg / (Angstrom s cm2)" + flux_unit_str = "erg / (Angstrom s cm2 pix2)" assert label_mouseover.as_text() == (f'Pixel x=00.0 y=00.0 Value +1.00000e+00 1e-17 {flux_unit_str}', # noqa 'World 13h41m46.5994s +26d59m58.6136s (ICRS)', '205.4441642302 26.9996148973 (deg)') @@ -128,7 +128,7 @@ def test_spectrum3d_parse(image_cube_hdu_obj, cubeviz_helper): label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 0, 'y': 0}}) - flux_unit_str = "erg / (Angstrom s cm2)" + flux_unit_str = "erg / (Angstrom s cm2 pix2)" assert label_mouseover.as_text() == (f'Pixel x=00.0 y=00.0 Value +1.00000e+00 1e-17 {flux_unit_str}', # noqa 'World 13h41m46.5994s +26d59m58.6136s (ICRS)', '205.4441642302 26.9996148973 (deg)') @@ -152,9 +152,10 @@ def test_spectrum3d_no_wcs_parse(cubeviz_helper): assert data.shape == (2, 3, 4) # x, y, z assert isinstance(data.coords, PaddedSpectrumWCS) assert_array_equal(flux.data, 1) - assert flux.units == 'nJy' + assert flux.units == 'nJy / pix2' +@pytest.mark.skip(reason="unskip after 3192 merged") def test_spectrum1d_parse(spectrum1d, cubeviz_helper): cubeviz_helper.load_data(spectrum1d) diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_tools.py b/jdaviz/configs/cubeviz/plugins/tests/test_tools.py index 8f08d55319..f5db618cd8 100644 --- a/jdaviz/configs/cubeviz/plugins/tests/test_tools.py +++ b/jdaviz/configs/cubeviz/plugins/tests/test_tools.py @@ -68,8 +68,22 @@ def test_spectrum_at_spaxel(cubeviz_helper, spectrum1d_cube_with_uncerts): assert uncert_viewer.toolbar.active_tool._mark.visible is True -def test_spectrum_at_spaxel_altkey_true(cubeviz_helper, spectrum1d_cube): - cubeviz_helper.load_data(spectrum1d_cube, data_label='test') +@pytest.mark.parametrize("cube_type", ["Surface Brightness", "Flux"]) +def test_spectrum_at_spaxel_altkey_true(cubeviz_helper, spectrum1d_cube, + spectrum1d_cube_sb_unit, cube_type): + + # test is parameterize to test a cube that is in Jy / sr (Surface Brightness) + # as well as Jy (Flux), to test that flux cubes, which are converted in the + # parser to flux / pix^2 surface brightness cubes, both work correctly. + + if cube_type == 'Surface Brightness': + cube = spectrum1d_cube_sb_unit + cube_unit = 'Jy / sr' + elif cube_type == 'Flux': + cube = spectrum1d_cube + cube_unit = 'Jy / pix2' + + cubeviz_helper.load_data(cube, data_label='test') flux_viewer = cubeviz_helper.app.get_viewer("flux-viewer") uncert_viewer = cubeviz_helper.app.get_viewer("uncert-viewer") @@ -88,7 +102,7 @@ def test_spectrum_at_spaxel_altkey_true(cubeviz_helper, spectrum1d_cube): label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 2, 'y': 1}}) - assert label_mouseover.as_text() == ('Pixel x=02.0 y=01.0 Value +1.40000e+01 Jy', + assert label_mouseover.as_text() == (f'Pixel x=02.0 y=01.0 Value +1.40000e+01 {cube_unit}', 'World 13h39m59.9192s +27d00m00.7200s (ICRS)', '204.9996633015 27.0001999996 (deg)') @@ -122,7 +136,7 @@ def test_spectrum_at_spaxel_altkey_true(cubeviz_helper, spectrum1d_cube): # Make sure coordinate info panel did not change label_mouseover._viewer_mouse_event(flux_viewer, {'event': 'mousemove', 'domain': {'x': 1, 'y': 1}}) - assert label_mouseover.as_text() == ('Pixel x=01.0 y=01.0 Value +1.30000e+01 Jy', + assert label_mouseover.as_text() == (f'Pixel x=01.0 y=01.0 Value +1.30000e+01 {cube_unit}', 'World 13h39m59.9461s +27d00m00.7200s (ICRS)', '204.9997755344 27.0001999998 (deg)') diff --git a/jdaviz/configs/default/plugins/markers/tests/test_markers_plugin.py b/jdaviz/configs/default/plugins/markers/tests/test_markers_plugin.py index 8d4e0af44e..d2c445b310 100644 --- a/jdaviz/configs/default/plugins/markers/tests/test_markers_plugin.py +++ b/jdaviz/configs/default/plugins/markers/tests/test_markers_plugin.py @@ -26,6 +26,8 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): cubeviz_helper.load_data(spectrum1d_cube, "test") fv = cubeviz_helper.app.get_viewer('flux-viewer') sv = cubeviz_helper.app.get_viewer('spectrum-viewer') + sb_unit = 'Jy / pix2' # cubes loaded in Jy have sb unit of Jy / pix2 + flux_unit = 'Jy' label_mouseover = cubeviz_helper.app.session.application._tools['g-coords-info'] @@ -41,7 +43,7 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): {'event': 'mousemove', 'domain': {'x': 0, 'y': 0}}) - assert label_mouseover.as_text() == ('Pixel x=00.0 y=00.0 Value +8.00000e+00 Jy', + assert label_mouseover.as_text() == (f'Pixel x=00.0 y=00.0 Value +8.00000e+00 {sb_unit}', 'World 13h39m59.9731s +27d00m00.3600s (ICRS)', '204.9998877673 27.0001000000 (deg)') @@ -60,7 +62,7 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): 'pixel_y': 0, 'pixel:unreliable': False, 'value': 8.0, - 'value:unit': 'Jy', + 'value:unit': sb_unit, 'value:unreliable': False}) mp._obj._on_viewer_key_event(fv, {'event': 'keydown', @@ -70,14 +72,14 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): # test event in spectrum viewer (with auto layer detection) # x = [4.62280007e-07, 4.62360028e-07] - # y = [28, 92] Jy + # y = [28, 92] Jy / pix2 label_mouseover._viewer_mouse_event(sv, {'event': 'mousemove', 'domain': {'x': 4.623e-7, 'y': 0}}) - assert label_mouseover.as_text() == ('Cursor 4.62300e-07, 0.00000e+00 Value +8.00000e+00 Jy', + assert label_mouseover.as_text() == (f'Cursor 4.62300e-07, 0.00000e+00 Value +8.00000e+00 {sb_unit}', # noqa 'Wave 4.62280e-07 m (0 pix)', - 'Flux 2.80000e+01 Jy') + f'Flux 2.80000e+01 {flux_unit}') assert label_mouseover.as_dict() == {'data_label': 'Spectrum (sum)', 'axes_x': 4.622800069238093e-07, 'axes_x:unit': 'm', @@ -85,9 +87,9 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): 'spectral_axis': 4.622800069238093e-07, 'spectral_axis:unit': 'm', 'axes_y': 28.0, - 'axes_y:unit': 'Jy', + 'axes_y:unit': flux_unit, 'value': 28.0, - 'value:unit': 'Jy'} + 'value:unit': flux_unit} mp._obj._on_viewer_key_event(sv, {'event': 'keydown', 'key': 'm'}) @@ -101,17 +103,17 @@ def test_markers_cubeviz(tmp_path, cubeviz_helper, spectrum1d_cube): {'event': 'mousemove', 'domain': {'x': 4.623e-7, 'y': 0}}) - assert label_mouseover.as_text() == ('Cursor 4.62300e-07, 0.00000e+00 Value +8.00000e+00 Jy', + assert label_mouseover.as_text() == (f'Cursor 4.62300e-07, 0.00000e+00 Value +8.00000e+00 {sb_unit}', # noqa '', '') assert label_mouseover.as_dict() == {'axes_x': 4.623e-07, 'axes_x:unit': 'm', 'axes_y': 0, - 'axes_y:unit': 'Jy', + 'axes_y:unit': flux_unit, 'data_label': '', 'spectral_axis': 4.623e-07, 'spectral_axis:unit': 'm', 'value': 0, - 'value:unit': 'Jy'} + 'value:unit': flux_unit} mp._obj._on_viewer_key_event(sv, {'event': 'keydown', 'key': 'm'}) diff --git a/jdaviz/configs/default/plugins/model_fitting/model_fitting.py b/jdaviz/configs/default/plugins/model_fitting/model_fitting.py index bdf882824c..3ed30d20f6 100644 --- a/jdaviz/configs/default/plugins/model_fitting/model_fitting.py +++ b/jdaviz/configs/default/plugins/model_fitting/model_fitting.py @@ -25,6 +25,7 @@ from jdaviz.configs.default.plugins.model_fitting.initializers import (MODELS, initialize, get_model_parameters) +from jdaviz.utils import _eqv_flux_to_sb_pixel __all__ = ['ModelFitting'] @@ -534,7 +535,10 @@ def _initialize_model_component(self, model_comp, comp_label, poly_order=None): init_x = masked_spectrum.spectral_axis init_y = masked_spectrum.flux - init_y = init_y.to(self._units['y'], u.spectral_density(init_x)) + # equivs for spectral density and flux<>flux/pix2. revisit + # when generalizing plugin UC equivs. + equivs = _eqv_flux_to_sb_pixel() + [u.spectral_density(init_x)] + init_y = init_y.to(self._units['y'], equivs) initialized_model = initialize( MODELS[model_comp](name=comp_label, diff --git a/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py b/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py index 41926d6330..9e701412ae 100644 --- a/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py +++ b/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py @@ -76,14 +76,23 @@ def test_model_ids(cubeviz_helper, spectral_cube_wcs): plugin.vue_add_model({}) -@pytest.mark.skip(reason="Needs #3156 after merging #3190") def test_parameter_retrieval(cubeviz_helper, spectral_cube_wcs): + + flux_unit = u.nJy + sb_unit = flux_unit / (u.pix * u.pix) + wav_unit = u.Hz + flux = np.ones((3, 4, 5)) flux[2, 2, :] = [1, 2, 3, 4, 5] - cubeviz_helper.load_data(Spectrum1D(flux=flux * u.nJy, wcs=spectral_cube_wcs), + cubeviz_helper.load_data(Spectrum1D(flux=flux * flux_unit, wcs=spectral_cube_wcs), data_label='test') + plugin = cubeviz_helper.plugins["Model Fitting"] + + # since cube_fit is true, output fit should be in SB units + # even though the spectral y axis is in 'flux' by default plugin.cube_fit = True + plugin.create_model_component("Linear1D", "L") with warnings.catch_warnings(): warnings.filterwarnings('ignore', message='Model is linear in parameters.*') @@ -92,14 +101,15 @@ def test_parameter_retrieval(cubeviz_helper, spectral_cube_wcs): params = cubeviz_helper.get_model_parameters() slope_res = np.zeros((3, 4)) slope_res[2, 2] = 1.0 - slope_res = slope_res * u.nJy / u.Hz + + slope_res = slope_res * sb_unit / wav_unit intercept_res = np.ones((3, 4)) intercept_res[2, 2] = 0 - intercept_res = intercept_res * u.nJy + intercept_res = intercept_res * sb_unit assert_quantity_allclose(params['model']['slope'], slope_res, - atol=1e-10 * u.nJy / u.Hz) + atol=1e-10 * sb_unit / wav_unit) assert_quantity_allclose(params['model']['intercept'], intercept_res, - atol=1e-10 * u.nJy) + atol=1e-10 * sb_unit) @pytest.mark.parametrize('unc', ('zeros', None)) @@ -225,6 +235,8 @@ def test_cube_fitting_backend(cubeviz_helper, unc, tmp_path): assert isinstance(fitted_spectrum, Spectrum1D) assert len(fitted_spectrum.shape) == 3 assert fitted_spectrum.shape == (IMAGE_SIZE_X, IMAGE_SIZE_Y, SPECTRUM_SIZE) + + # since were not loading data into app, results are just u.Jy not u.Jy / pix2 assert fitted_spectrum.flux.unit == u.Jy assert not np.all(fitted_spectrum.flux.value == 0) @@ -275,7 +287,9 @@ def test_cube_fitting_backend(cubeviz_helper, unc, tmp_path): data_sci = cubeviz_helper.app.data_collection["fitted_cube.fits[SCI]"] flux_sci = data_sci.get_component("flux") assert_allclose(flux_sci.data, fitted_spectrum.flux.value) - assert flux_sci.units == flux_unit_str + # now that the flux cube was loaded into cubeviz, there will be a factor + # of pix2 applied to the flux unit + assert flux_sci.units == f'{flux_unit_str} / pix2' coo = data_sci.coords.pixel_to_world(1, 0, 2) assert_allclose(coo[0].value, coo_expected[0].value) # SpectralCoord assert_allclose([coo[1].ra.deg, coo[1].dec.deg], @@ -382,7 +396,6 @@ def test_incompatible_units(specviz_helper, spectrum1d): mf.calculate_fit(add_data=True) -@pytest.mark.skip(reason="Needs #3156 after merging #3190") def test_cube_fit_with_nans(cubeviz_helper): flux = np.ones((7, 8, 9)) * u.nJy flux[:, :, 0] = np.nan @@ -403,7 +416,6 @@ def test_cube_fit_with_nans(cubeviz_helper): assert mf._obj.component_models[0]['compat_display_units'] is False -@pytest.mark.skip(reason="Needs #3156 after merging #3190") def test_cube_fit_with_subset_and_nans(cubeviz_helper): # Also test with existing mask flux = np.ones((7, 8, 9)) * u.nJy diff --git a/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py b/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py index b75551cbc8..a5bc92f011 100644 --- a/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py +++ b/jdaviz/configs/default/plugins/model_fitting/tests/test_plugin.py @@ -134,7 +134,6 @@ def test_register_model_uncertainty_is_none(specviz_helper, spectrum1d): assert np.allclose(param["std"], expected_uncertainties[param["name"]], rtol=0.01) -@pytest.mark.skip(reason="Needs #3156 after merging #3190") def test_register_cube_model(cubeviz_helper, spectrum1d_cube): with warnings.catch_warnings(): warnings.simplefilter('ignore') @@ -148,6 +147,7 @@ def test_register_cube_model(cubeviz_helper, spectrum1d_cube): # changing the lable should set auto to False, but the event may not have triggered yet modelfit_plugin._obj.results_label_auto = False modelfit_plugin.cube_fit = True + modelfit_plugin.reestimate_model_parameters() assert modelfit_plugin._obj.results_label_default == 'model' assert modelfit_plugin._obj.results_label == test_label with warnings.catch_warnings(): @@ -156,7 +156,6 @@ def test_register_cube_model(cubeviz_helper, spectrum1d_cube): assert test_label in cubeviz_helper.app.data_collection -@pytest.mark.skip(reason="Needs #3156 after merging #3190") def test_fit_cube_no_wcs(cubeviz_helper): # This is like when user do something to a cube outside of Jdaviz # and then load it back into a new instance of Cubeviz for further analysis. diff --git a/jdaviz/configs/default/plugins/viewers.py b/jdaviz/configs/default/plugins/viewers.py index 5bf05593e6..ff9077542d 100644 --- a/jdaviz/configs/default/plugins/viewers.py +++ b/jdaviz/configs/default/plugins/viewers.py @@ -31,7 +31,9 @@ from jdaviz.core.template_mixin import WithCache from jdaviz.core.user_api import ViewerUserApi from jdaviz.utils import (ColorCycler, get_subset_type, _wcs_only_label, - layer_is_image_data, layer_is_not_dq) + layer_is_image_data, layer_is_not_dq, + _eqv_sb_per_pixel_to_per_angle) +from jdaviz.core.validunits import check_if_unit_is_per_solid_angle uc = UnitConverter() @@ -754,24 +756,43 @@ def set_plot_axes(self): u.bol, u.AB, u.ST ] - locally_defined_sb_units = [ - unit / u.sr for unit in locally_defined_flux_units - ] - - if any(y_unit.is_equivalent(unit) for unit in locally_defined_sb_units): - flux_unit_type = "Surface Brightness" - elif any(y_unit.is_equivalent(unit) for unit in locally_defined_flux_units): - flux_unit_type = 'Flux' - elif ( - y_unit.is_equivalent(u.DN) or - y_unit.is_equivalent(u.electron / u.s) or - y_unit.physical_type == 'dimensionless' - ): - # electron / s or 'dimensionless_unscaled' should be labeled counts - flux_unit_type = "Counts" - elif y_unit.is_equivalent(u.W): - flux_unit_type = "Luminosity" + # get square angle from 'sb' display unit + sb_unit = self.jdaviz_app._get_display_unit(axis='sb') + if sb_unit is not None: + solid_angle_unit = check_if_unit_is_per_solid_angle(sb_unit, return_unit=True) else: + solid_angle_unit = None + + # if solid angle is present in denominator, check physical type of numerator + # if numerator is a flux type the display unit is a 'surface brightness', otherwise + # default to the catchall 'flux density' label + flux_unit_type = None + + if solid_angle_unit is not None: + + for un in locally_defined_flux_units: + locally_defined_sb_unit = un / solid_angle_unit + + # create an equivalency for each flux unit for flux <> flux/pix2. + # for similar reasons to the 'untranslatable units' issue, custom + # equivs. can't be combined, so a workaround is creating an eqiv + # for each flux that may need an additional equiv. + angle_to_pixel_equiv = _eqv_sb_per_pixel_to_per_angle(un) + + if y_unit.is_equivalent(locally_defined_sb_unit, angle_to_pixel_equiv): + flux_unit_type = "Surface Brightness" + elif y_unit.is_equivalent(un): + flux_unit_type = 'Flux' + elif y_unit.is_equivalent(u.electron / u.s) or y_unit.physical_type == 'dimensionless': # noqa + # electron / s or 'dimensionless_unscaled' should be labeled counts + flux_unit_type = "Counts" + elif y_unit.is_equivalent(u.W): + flux_unit_type = "Luminosity" + if flux_unit_type is not None: + # if we determined a label, stop checking + break + + if flux_unit_type is None: # default to Flux Density for flux density or uncaught types flux_unit_type = "Flux density" diff --git a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py index 97e3f90dd1..db21424a24 100644 --- a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py +++ b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py @@ -27,6 +27,8 @@ __all__ = ['SimpleAperturePhotometry'] +PIX2 = u.pix * u.pix # define square pixel unit which is used around the plugin + @tray_registry('imviz-aper-phot-simple', label="Aperture Photometry") class SimpleAperturePhotometry(PluginTemplateMixin, ApertureSubsetSelectMixin, @@ -72,7 +74,13 @@ class SimpleAperturePhotometry(PluginTemplateMixin, ApertureSubsetSelectMixin, # surface brightness display unit display_unit = Unicode("").tag(sync=True) - # flux scaling display unit will always be flux, not sb + # angle component of `display_unit`, to avoid repetition of seperating + # it out from `display_unit` + + display_solid_angle_unit = Unicode("").tag(sync=True) + + # flux scaling display unit will always be flux, not sb. again its own + # traitlet to avoid avoid repetition of seperating it out from `display_unit` flux_scaling_display_unit = Unicode("").tag(sync=True) def __init__(self, *args, **kwargs): @@ -112,16 +120,23 @@ def __init__(self, *args, **kwargs): # Custom dataset filters for Cubeviz if self.config == "cubeviz": def valid_cubeviz_datasets(data): + comp = data.get_component(data.main_components[0]) img_unit = u.Unit(comp.units) if comp.units else u.dimensionless_unscaled + solid_angle_unit = check_if_unit_is_per_solid_angle(img_unit, return_unit=True) + if solid_angle_unit is None: # this is encountered sometimes ?? + return + + # multiply out solid angle so we can check physical type of numerator + img_unit *= solid_angle_unit + acceptable_types = ['spectral flux density wav', 'photon flux density wav', 'spectral flux density', - 'photon flux density', - 'surface brightness'] + 'photon flux density'] + return ((data.ndim in (2, 3)) and - ((img_unit == (u.MJy / u.sr)) or - (img_unit.physical_type in acceptable_types))) + (img_unit.physical_type in acceptable_types)) self.dataset.add_filter(valid_cubeviz_datasets) self.session.hub.subscribe(self, SliceValueUpdatedMessage, @@ -195,6 +210,9 @@ def _on_display_units_changed(self, event={}): new_unit = u.Unit(self.display_unit) bg = self.background_value * prev_unit + + # will need to add additonal custom equiv here once + # pix2<>angle is enabled self.background_value = bg.to_value( new_unit, u.spectral_density(self._cube_wave)) @@ -216,36 +234,41 @@ def _set_display_unit_of_selected_dataset(self): surface brightness. """ - if not self.dataset_selected or not self.aperture_selected: + # all cubes are in sb so we can get display unit for plugin from SB display unit + # this can be changed to listen specifically to changes in surface brightness + # from UC plugin GlobalDisplayUnitChange message, but wiill require some refactoring + disp_unit = self.app._get_display_unit('sb') + + # this check needs to be here because 'get_display_unit' will sometimes + # return non surface brightness units or even None when the app is starting + # up. this can be removed once that is fixed (see PR #3144) + if disp_unit is None or not check_if_unit_is_per_solid_angle(disp_unit): self.display_unit = '' self.flux_scaling_display_unit = '' return - data = self.dataset.selected_dc_item - if isinstance(data, list): - data = data[0] - comp = data.get_component(data.main_components[0]) - if comp.units: - # if data is something-per-solid-angle, its a SB unit and we should - # use the selected global display unit for SB - if check_if_unit_is_per_solid_angle(comp.units): - display_unit_type = 'sb' - else: - display_unit_type = 'flux' + self.display_unit = disp_unit - disp_unit = self.app._get_display_unit(display_unit_type) + # get angle componant of surface brightness + # note: could add 'axis=angle' when cleaning this code up to avoid repeating this - self.display_unit = disp_unit + display_solid_angle_unit = check_if_unit_is_per_solid_angle(disp_unit, return_unit=True) + if display_solid_angle_unit is not None: + self.display_solid_angle_unit = display_solid_angle_unit.to_string() + else: + # there should always be a solid angle, but i think this is + # encountered sometimes when initializing something.. + self.display_solid_angle_unit = '' - # now get display unit for flux_scaling_display_unit. this unit will always - # be in flux, but it will not be derived from the global flux display unit - # note : need to generalize this for non-sr units eventually - fs_unit = u.Unit(disp_unit) * u.sr - self.flux_scaling_display_unit = fs_unit.to_string() + # flux scaling will be applied when the solid angle componant is + # multiplied out, so use 'flux' display unit + fs_unit = self.app._get_display_unit('flux') + self.flux_scaling_display_unit = fs_unit - else: - self.display_unit = '' - self.flux_scaling_display_unit = '' + # if cube loaded is per-pixel-squared sb (i.e flux cube loaded) + # pixel_area should be fixed to 1 + if self.display_solid_angle_unit == 'pix2': + self.pixel_area = 1.0 def _get_defaults_from_metadata(self, dataset=None): defaults = {} @@ -353,7 +376,7 @@ def _dataset_selected_changed(self, event={}): # get correct display unit for newly selected dataset if self.config == 'cubeviz': - # set display_unit and flux_scaling_display_unit + # sets display_unit and flux_scaling_display_unit traitlets self._set_display_unit_of_selected_dataset() # auto-populate background, if applicable. @@ -505,7 +528,7 @@ def calculate_photometry(self, dataset=None, aperture=None, background=None, Background to subtract, same unit as data. Automatically computed if ``background`` is set to a subset. pixel_area : float, optional - Pixel area in arcsec squared, only used if sr in data unit. + Pixel area in arcsec squared, only used if data unit is a surface brightness unit. counts_factor : float, optional Factor to convert data unit to counts, in unit of flux/counts. flux_scaling : float, optional @@ -624,7 +647,8 @@ def calculate_photometry(self, dataset=None, aperture=None, background=None, ycenter = reg.center.y if data.coords is not None: if self.config == "cubeviz" and data.ndim > 2: - sky_center = w.pixel_to_world(self._cubeviz_slice_ind, ycenter, xcenter)[1] + sky_center = w.pixel_to_world(self._cubeviz_slice_ind, + ycenter, xcenter)[1] else: # "imviz" sky_center = w.pixel_to_world(xcenter, ycenter) else: @@ -667,9 +691,11 @@ def calculate_photometry(self, dataset=None, aperture=None, background=None, # data units and does not need to be converted if ((self.config == 'cubeviz') and (flux_scaling is None) and (self.flux_scaling is not None)): - # update eventaully to handle non-sr SB units + + # convert flux_scaling from flux display unit to native flux unit flux_scaling = (self.flux_scaling * u.Unit(self.flux_scaling_display_unit)).to_value( # noqa: E501 - img_unit * u.sr, u.spectral_density(self._cube_wave)) + img_unit * self.display_solid_angle_unit, + u.spectral_density(self._cube_wave)) try: flux_scale = float(flux_scaling if flux_scaling is not None else self.flux_scaling) @@ -693,11 +719,25 @@ def calculate_photometry(self, dataset=None, aperture=None, background=None, rawsum = phot_table['sum'][0] if include_pixarea_fac: - pixarea = pixarea * (u.arcsec * u.arcsec / (u.pix * u.pix)) - # NOTE: Sum already has npix value encoded, so we simply apply the npix unit here. + # convert pixarea, which is in arcsec2/pix2 to the display solid angle unit / pix2 + + if self.config == 'imviz': + # can remove once unit conversion implemented in imviz and + # display_solid_angle_unit traitlet is set, for now it will always be the data units + display_solid_angle_unit = check_if_unit_is_per_solid_angle(comp.units, + return_unit=True) + + elif self.config == 'cubeviz': + display_solid_angle_unit = u.Unit(self.display_solid_angle_unit) + + # if angle unit is pix2, pixarea should be 1 pixel2 per pixel2 + if display_solid_angle_unit == PIX2: + pixarea_fac = 1 * PIX2 + else: + pixarea = pixarea * (u.arcsec * u.arcsec / PIX2) + # NOTE: Sum already has npix value encoded, so we simply apply the npix unit here. + pixarea_fac = PIX2 * pixarea.to(display_solid_angle_unit / PIX2) - # note: need to generalize this to non-steradian surface brightness units - pixarea_fac = (u.pix * u.pix) * pixarea.to(u.sr / (u.pix * u.pix)) phot_table['sum'] = [rawsum * pixarea_fac] else: pixarea_fac = None @@ -1227,7 +1267,8 @@ def _curve_of_growth(data, centroid, aperture, final_sum, wcs=None, background=0 else: sum_unit = None if sum_unit and pixarea_fac is not None: - sum_unit *= pixarea_fac.unit + # multiply data unit by its solid angle to convert sum in sb to sum in flux + sum_unit *= check_if_unit_is_per_solid_angle(sum_unit, return_unit=True) if hasattr(aperture, 'to_pixel'): aperture = aperture.to_pixel(wcs) diff --git a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue index ba02f49305..1041ab1768 100644 --- a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue +++ b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue @@ -95,7 +95,8 @@ persistent-hint /> - + + surface brightness. + """ + cube = spectrum1d_cube_custom_fluxunit(fluxunit=u.Jy / sq_angle_unit, + shape=(25, 25, 4), with_uncerts=True) + cubeviz_helper.load_data(cube, data_label="Test Cube") + + uc = cubeviz_helper.plugins['Unit Conversion'] + assert uc.spectral_y_type == 'Flux' # initial selection should be Flux + + plugin = cubeviz_helper.app.get_tray_item_from_name('specviz-line-analysis') + plugin.keep_active = True + + plugin.dataset_selected = 'Spectrum (sum)' + plugin.continuum_subset_selected = 'Surrounding' + + results = plugin.results + assert results[0]['unit'] == 'W / m2' + + # now change to surface brightness + uc.spectral_y_type = 'Surface Brightness' + + results = plugin.results + line_flux_before_unit_conversion = results[0] + # convert back and forth between unit<>str for string format consistency + unit_str = u.Unit(f'W / {sq_angle_unit.to_string()} m2').to_string() + assert line_flux_before_unit_conversion['unit'] == unit_str + + # change flux unit and make sure result stays the same after conversion + uc.flux_unit.selected = 'MJy' + results = plugin.results + np.testing.assert_allclose(float(results[0]['result']), + float(line_flux_before_unit_conversion['result'])) + np.testing.assert_allclose(float(results[0]['uncertainty']), + float(line_flux_before_unit_conversion['uncertainty'])) + + def test_user_api(specviz_helper, spectrum1d): label = "Test 1D Spectrum" specviz_helper.load_data(spectrum1d, data_label=label) diff --git a/jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py b/jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py index f6607d745d..b967e90ffe 100644 --- a/jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py +++ b/jdaviz/configs/specviz/plugins/unit_conversion/tests/test_unit_conversion.py @@ -137,16 +137,23 @@ def test_non_stddev_uncertainty(specviz_helper): ) -def test_unit_translation(cubeviz_helper): - # custom cube so PIXAR_SR is in metadata, and Flux units, and in MJy +def cubeviz_wcs_dict(): + # returns a WCS obj and dictionary used for cubeviz tests here wcs_dict = {"CTYPE1": "WAVE-LOG", "CTYPE2": "DEC--TAN", "CTYPE3": "RA---TAN", "CRVAL1": 4.622e-7, "CRVAL2": 27, "CRVAL3": 205, "CDELT1": 8e-11, "CDELT2": 0.0001, "CDELT3": -0.0001, "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0, "PIXAR_SR": 8e-11} w = WCS(wcs_dict) + return w, wcs_dict + + +@pytest.mark.parametrize("angle_unit", [u.sr, u.pix*u.pix]) +def test_unit_translation(cubeviz_helper, angle_unit): + # custom cube so PIXAR_SR is in metadata, and Flux units, and in MJy + 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.MJy / u.sr, wcs=w, meta=wcs_dict) + cube = Spectrum1D(flux=flux * u.MJy / angle_unit, wcs=w, meta=wcs_dict) cubeviz_helper.load_data(cube, data_label="test") center = PixCoord(5, 10) @@ -176,23 +183,23 @@ def test_unit_translation(cubeviz_helper): y_display_unit = u.Unit(viewer_1d.state.y_display_unit) # check if units translated - assert y_display_unit == u.MJy / u.sr + assert y_display_unit == u.MJy / angle_unit # get_data(use_display_units=True) should return surface brightness-like units - assert cubeviz_helper.app._get_display_unit('spectral_y') == u.MJy / u.sr - assert cubeviz_helper.get_data('Spectrum (sum)', use_display_units=True).unit == u.MJy / u.sr + assert cubeviz_helper.app._get_display_unit('spectral_y') == u.MJy / angle_unit + assert cubeviz_helper.get_data('Spectrum (sum)', use_display_units=True).unit == u.MJy / angle_unit # noqa + + +@pytest.mark.parametrize("angle_unit", [u.sr, u.pix*u.pix]) +def test_sb_unit_conversion(cubeviz_helper, angle_unit): + angle_str = angle_unit.to_string() -def test_sb_unit_conversion(cubeviz_helper): # custom cube to have Surface Brightness units - wcs_dict = {"CTYPE1": "WAVE-LOG", "CTYPE2": "DEC--TAN", "CTYPE3": "RA---TAN", - "CRVAL1": 4.622e-7, "CRVAL2": 27, "CRVAL3": 205, - "CDELT1": 8e-11, "CDELT2": 0.0001, "CDELT3": -0.0001, - "CRPIX1": 0, "CRPIX2": 0, "CRPIX3": 0, "PIXAR_SR": 8e-11} - w = WCS(wcs_dict) + 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.MJy / u.sr), wcs=w, meta=wcs_dict) + cube = Spectrum1D(flux=flux * (u.MJy / angle_unit), wcs=w, meta=wcs_dict) cubeviz_helper.load_data(cube, data_label="test") uc_plg = cubeviz_helper.plugins['Unit Conversion'] @@ -202,6 +209,9 @@ def test_sb_unit_conversion(cubeviz_helper): assert uc_plg.spectral_y_type == 'Flux' # flux choices is populated with flux units assert uc_plg.flux_unit.choices + # and angle choices should be the only the input angle + assert len(uc_plg.angle_unit.choices) == 1 + assert angle_str in uc_plg.angle_unit.choices # to have access to display units viewer_1d = cubeviz_helper.app.get_viewer( @@ -212,7 +222,7 @@ def test_sb_unit_conversion(cubeviz_helper): # Surface Brightness conversion uc_plg.flux_unit = 'Jy' y_display_unit = u.Unit(viewer_1d.state.y_display_unit) - assert y_display_unit == u.Jy / u.sr + assert y_display_unit == u.Jy / angle_unit label_mouseover = cubeviz_helper.app.session.application._tools["g-coords-info"] flux_viewer = cubeviz_helper.app.get_viewer( cubeviz_helper._default_flux_viewer_reference_name @@ -221,21 +231,28 @@ def test_sb_unit_conversion(cubeviz_helper): flux_viewer, {"event": "mousemove", "domain": {"x": 10, "y": 8}} ) assert label_mouseover.as_text() == ( - "Pixel x=00010.0 y=00008.0 Value +1.00000e+06 Jy / sr", + f"Pixel x=00010.0 y=00008.0 Value +1.00000e+06 Jy / {angle_str}", "World 13h39m59.7037s +27d00m03.2400s (ICRS)", "204.9987654313 27.0008999946 (deg)") # Try a second conversion uc_plg.flux_unit = 'W / Hz m2' + + if angle_unit == u.pix * u.pix: # unit string order is different for pix2 vs sr + str_unit = 'W / (Hz m2 pix2)' + elif angle_unit == u.sr: + str_unit = 'W / (Hz sr m2)' + y_display_unit = u.Unit(viewer_1d.state.y_display_unit) - assert y_display_unit == u.Unit("W / (Hz sr m2)") + assert y_display_unit == u.Unit(str_unit) y_display_unit = u.Unit(viewer_1d.state.y_display_unit) label_mouseover._viewer_mouse_event( flux_viewer, {"event": "mousemove", "domain": {"x": 10, "y": 8}} ) + assert label_mouseover.as_text() == ( - "Pixel x=00010.0 y=00008.0 Value +1.00000e-20 W / (Hz sr m2)", + f"Pixel x=00010.0 y=00008.0 Value +1.00000e-20 {str_unit}", "World 13h39m59.7037s +27d00m03.2400s (ICRS)", "204.9987654313 27.0008999946 (deg)") @@ -248,7 +265,7 @@ def test_sb_unit_conversion(cubeviz_helper): flux_viewer, {"event": "mousemove", "domain": {"x": 10, "y": 8}} ) assert label_mouseover.as_text() == ( - "Pixel x=00010.0 y=00008.0 Value +1.00000e+00 MJy / sr", + f"Pixel x=00010.0 y=00008.0 Value +1.00000e+00 MJy / {angle_str}", "World 13h39m59.7037s +27d00m03.2400s (ICRS)", "204.9987654313 27.0008999946 (deg)") @@ -279,3 +296,58 @@ def test_contour_unit_conversion(cubeviz_helper, spectrum1d_cube_fluxunit_jy_per uc_plg.flux_unit = 'MJy' assert np.allclose(po_plg.contour_max.value, 1.99e-4) + + +@pytest.mark.parametrize("angle_unit", [u.sr, u.pix*u.pix]) +def test_cubeviz_flux_sb_translation_counts(cubeviz_helper, angle_unit): + + """ + When a cube is loaded in counts, 'count' should be the only + available option for flux unit. The y axis can be translated + between flux and sb. Test a flux cube which will be converted + to ct/pix2, and a sb cube ct/sr. + """ + + angle_str = angle_unit.to_string() + + # custom cube to have Surface Brightness 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.ct / angle_unit), wcs=w, meta=wcs_dict) + cubeviz_helper.load_data(cube, data_label="test") + + uc_plg = cubeviz_helper.plugins['Unit Conversion'] + uc_plg.open_in_tray() + + # ensure that per solid angle cube defaults to Flux spectrum + assert uc_plg.spectral_y_type == 'Flux' + # flux choices is populated with only one choice, counts + assert len(uc_plg.flux_unit.choices) == 1 + assert 'ct' in uc_plg.flux_unit.choices + # and angle choices should be the only the input angle + assert len(uc_plg.angle_unit.choices) == 1 + assert angle_str in uc_plg.angle_unit.choices + + # to have access to display units + viewer_1d = cubeviz_helper.app.get_viewer( + cubeviz_helper._default_spectrum_viewer_reference_name) + + # do a spectral y axis translation from Flux to Surface Brightness + uc_plg.spectral_y_type.selected = 'Surface Brightness' + + y_display_unit = u.Unit(viewer_1d.state.y_display_unit) + assert y_display_unit == u.ct / angle_unit + + # and test mouseover info + label_mouseover = cubeviz_helper.app.session.application._tools["g-coords-info"] + flux_viewer = cubeviz_helper.app.get_viewer( + cubeviz_helper._default_flux_viewer_reference_name + ) + label_mouseover._viewer_mouse_event( + flux_viewer, {"event": "mousemove", "domain": {"x": 10, "y": 8}} + ) + assert label_mouseover.as_text() == ( + 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)") diff --git a/jdaviz/configs/specviz/plugins/unit_conversion/unit_conversion.py b/jdaviz/configs/specviz/plugins/unit_conversion/unit_conversion.py index bca8ff5bde..ae7ec6093b 100644 --- a/jdaviz/configs/specviz/plugins/unit_conversion/unit_conversion.py +++ b/jdaviz/configs/specviz/plugins/unit_conversion/unit_conversion.py @@ -189,6 +189,7 @@ def _on_glue_y_display_unit_changed(self, y_unit_str): x_unit = u.Unit(self.spectral_unit.selected) y_unit_str = _valid_glue_display_unit(y_unit_str, self.spectrum_viewer, 'y') y_unit = u.Unit(y_unit_str) + y_unit_solid_angle = check_if_unit_is_per_solid_angle(y_unit_str, return_unit=True) if not check_if_unit_is_per_solid_angle(y_unit_str) and y_unit_str != self.flux_unit.selected: # noqa flux_choices = create_flux_equivalencies_list(y_unit, x_unit) @@ -201,10 +202,11 @@ def _on_glue_y_display_unit_changed(self, y_unit_str): # if the y-axis is set to surface brightness, # untranslatable units need to be removed from the flux choices - if check_if_unit_is_per_solid_angle(y_unit_str): - flux_choices = create_flux_equivalencies_list(y_unit * u.sr, x_unit) + if y_unit_solid_angle: + flux_choices = [(y_unit * y_unit_solid_angle).to_string()] + flux_choices += create_flux_equivalencies_list(y_unit * y_unit_solid_angle, x_unit) self.flux_unit.choices = flux_choices - flux_unit = str(y_unit * u.sr) + flux_unit = str(y_unit * y_unit_solid_angle) # We need to set the angle_unit before triggering _on_flux_unit_changed via # setting self.flux_unit.selected below, or the lack of angle unit will make it think # we're in Flux units. @@ -216,7 +218,11 @@ def _on_glue_y_display_unit_changed(self, y_unit_str): # sets the angle unit drop down and the surface brightness read-only text if self.app.data_collection[0]: dc_unit = self.app.data_collection[0].get_component("flux").units - self.angle_unit.choices = create_angle_equivalencies_list(dc_unit) + + # angle choices will be angle equivalencies to the solid-angle component of the cube + dc_solid_angle_unit = check_if_unit_is_per_solid_angle(dc_unit, return_unit=True) + + self.angle_unit.choices = create_angle_equivalencies_list(dc_solid_angle_unit) self.angle_unit.selected = self.angle_unit.choices[0] self.sb_unit_selected = self._append_angle_correctly( self.flux_unit.selected, @@ -238,7 +244,8 @@ def _on_glue_y_display_unit_changed(self, y_unit_str): if not self.flux_unit.selected: y_display_unit = self.spectrum_viewer.state.y_display_unit - self.flux_unit.selected = (str(u.Unit(y_display_unit * u.sr))) + flux_unit_str = str(u.Unit(y_display_unit * y_unit_solid_angle)) + self.flux_unit.selected = flux_unit_str @observe('spectral_unit_selected') def _on_spectral_unit_changed(self, *args): @@ -350,6 +357,7 @@ def _find_and_convert_contour_units(self, msg=None, yunit=None): layer.attribute_display_unit = yunit def _translate(self, spectral_y_type=None): + # currently unsupported, can be supported with a scale factor if self.app.config == 'specviz': return @@ -364,18 +372,21 @@ def _translate(self, spectral_y_type=None): if not self.flux_unit.choices: return + selected_display_solid_angle_unit = u.Unit(self.angle_unit_selected) + spec_axis_ang_unit = check_if_unit_is_per_solid_angle(spec_units) + # Surface Brightness -> Flux - if check_if_unit_is_per_solid_angle(spec_units) and spectral_y_type == 'Flux': - spec_units *= u.sr + if spec_axis_ang_unit and spectral_y_type == 'Flux': + spec_units *= selected_display_solid_angle_unit # update display units self.spectrum_viewer.state.y_display_unit = str(spec_units) # Flux -> Surface Brightness - elif (not check_if_unit_is_per_solid_angle(spec_units) - and spectral_y_type == 'Surface Brightness'): - spec_units /= u.sr + elif (not spec_axis_ang_unit and spectral_y_type == 'Surface Brightness'): + spec_units /= selected_display_solid_angle_unit # update display units self.spectrum_viewer.state.y_display_unit = str(spec_units) + # entered the translator when we shouldn't translate else: return @@ -388,7 +399,7 @@ def _translate(self, spectral_y_type=None): self.spectrum_viewer.reset_limits() def _append_angle_correctly(self, flux_unit, angle_unit): - if angle_unit not in ['pix', 'sr']: + if angle_unit not in ['pix2', 'sr']: self.sb_unit_selected = flux_unit return flux_unit if '(' in flux_unit: @@ -399,6 +410,9 @@ def _append_angle_correctly(self, flux_unit, angle_unit): sb_unit_selected = flux_unit + ' / ' + angle_unit if sb_unit_selected: - self.sb_unit_selected = sb_unit_selected + # convert string to and from u.Unit to get rid of any + # formatting inconstancies, order of units in string + # for a composite unit matters + sb_unit_selected = u.Unit(sb_unit_selected).to_string() return sb_unit_selected diff --git a/jdaviz/configs/specviz/tests/test_viewers.py b/jdaviz/configs/specviz/tests/test_viewers.py index ba706c32e1..cdbd6ea5ee 100644 --- a/jdaviz/configs/specviz/tests/test_viewers.py +++ b/jdaviz/configs/specviz/tests/test_viewers.py @@ -4,6 +4,7 @@ from specutils import Spectrum1D +@pytest.mark.skip(reason="unskip after 3192 merged") @pytest.mark.parametrize( ('input_unit', 'y_axis_label'), [(u.MJy, 'Flux'), diff --git a/jdaviz/conftest.py b/jdaviz/conftest.py index 0be26cb176..a5f2e11b8c 100644 --- a/jdaviz/conftest.py +++ b/jdaviz/conftest.py @@ -306,6 +306,13 @@ def spectrum1d_cube_fluxunit_jy_per_steradian(): with_uncerts=True) +@pytest.fixture +def spectrum1d_cube_sb_unit(): + # similar fixture to spectrum1d_cube_fluxunit_jy_per_steradian, but no uncerts + # and different shape. can probably remove one of these eventually + return _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy / u.sr) + + @pytest.fixture def mos_spectrum1d(mos_spectrum2d): ''' diff --git a/jdaviz/core/marks.py b/jdaviz/core/marks.py index 3585369a90..4bcd52d2e4 100644 --- a/jdaviz/core/marks.py +++ b/jdaviz/core/marks.py @@ -5,7 +5,7 @@ from bqplot.marks import Lines, Label, Scatter from glue.core import HubListener from specutils import Spectrum1D -from jdaviz.utils import _eqv_pixar_sr +from jdaviz.utils import _eqv_pixar_sr, _eqv_flux_to_sb_pixel from jdaviz.core.events import GlobalDisplayUnitChanged from jdaviz.core.events import (SliceToolStateMessage, LineIdentifyMessage, @@ -113,9 +113,14 @@ def set_y_unit(self, unit=None): if self.viewer.default_class is Spectrum1D: spec = self.viewer.state.reference_data.get_object(cls=Spectrum1D) eqv = u.spectral_density(spec.spectral_axis) + if ('_pixel_scale_factor' in spec.meta): eqv += _eqv_pixar_sr(spec.meta['_pixel_scale_factor']) - y = (self.y * self.yunit).to_value(unit, equivalencies=eqv) + + # add equiv for flux <> flux/pix2 + eqv += _eqv_flux_to_sb_pixel() + + y = (self.y * self.yunit).to_value(unit, equivalencies=eqv) else: y = (self.y * self.yunit).to_value(unit) self.yunit = unit diff --git a/jdaviz/core/validunits.py b/jdaviz/core/validunits.py index 578c94f175..cfafb0b08a 100644 --- a/jdaviz/core/validunits.py +++ b/jdaviz/core/validunits.py @@ -1,10 +1,14 @@ from astropy import units as u import itertools -__all__ = ['units_to_strings', 'create_spectral_equivalencies_list', +__all__ = ['supported_sq_angle_units', 'units_to_strings', 'create_spectral_equivalencies_list', 'create_flux_equivalencies_list', 'check_if_unit_is_per_solid_angle'] +def supported_sq_angle_units(): + return [u.pix*u.pix, u.sr] + + def units_to_strings(unit_list): """Convert equivalencies into readable versions of the units. @@ -73,9 +77,12 @@ def create_flux_equivalencies_list(flux_unit, spectral_axis_unit): # Get local flux units. locally_defined_flux_units = ['Jy', 'mJy', 'uJy', 'MJy', - 'W / (Hz m2)', 'eV / (Hz s m2)', - 'erg / (s cm2 Angstrom)', 'erg / (Hz s cm2)', - 'ph / (Angstrom s cm2)', 'ph / (Hz s cm2)' + 'W / (Hz m2)', + 'eV / (s m2 Hz)', + 'erg / (s cm2 Hz)', + 'erg / (s cm2 Angstrom)', + 'ph / (Angstrom s cm2)', + 'ph / (Hz s cm2)', ] local_units = [u.Unit(unit) for unit in locally_defined_flux_units] @@ -90,39 +97,63 @@ def create_flux_equivalencies_list(flux_unit, spectral_axis_unit): return sorted(units_to_strings(local_units)) + flux_unit_equivalencies_titles -def create_angle_equivalencies_list(unit): - # first, convert string to u.Unit obj. - # this will take care of some formatting consistency like - # turning something like Jy / (degree*degree) to Jy / deg**2 - # and erg sr^1 to erg / sr - if isinstance(unit, u.core.Unit) or isinstance(unit, u.core.CompositeUnit): - unit_str = unit.to_string() - elif isinstance(unit, str): - unit = u.Unit(unit) - unit_str = unit.to_string() - elif unit == 'ct': - return ['pix'] - else: - raise ValueError('Unit must be u.Unit, or string that can be converted into a u.Unit') +def create_angle_equivalencies_list(solid_angle_unit): - if '/' in unit_str: - # might be comprised of several units in denom. - denom = unit_str.split('/')[-1].split() - return denom - else: - # this could be where force / u.pix - return ['pix'] + """ + Return valid angles that `solid_angle_unit` (which should be a solid angle + physical type, or square pixel), can be translated to in the unit conversion + plugin. These options will populate the dropdown menu for 'angle unit' in + the Unit Conversion plugin. + + Parameters + ---------- + solid_angle_unit : str or u.Unit + Unit object or string representation of unit that is a 'solid angle' + or square pixel physical type. + + Returns + ------- + equivalent_angle_units : list of str + String representation of units that `solid_angle_unit` can be + translated to. + """ + + if solid_angle_unit is None: + # if there was no solid angle in the unit when calling this function + # can only represent that unit as per square pixel + return ['pix^2'] + + # cast to unit then back to string to account for formatting inconsistencies + # in strings that represent units + if isinstance(solid_angle_unit, str): + solid_angle_unit = u.Unit(solid_angle_unit) + unit_str = solid_angle_unit.to_string() + + # uncomment and expand this list once translating between solid + # angles and between solid angle and solid pixel is enabled + # equivalent_angle_units = ['sr', 'pix^2'] + equivalent_angle_units = [] + if unit_str not in equivalent_angle_units: + equivalent_angle_units += [unit_str] -def check_if_unit_is_per_solid_angle(unit): + return equivalent_angle_units + + +def check_if_unit_is_per_solid_angle(unit, return_unit=False): """ - Check if a given u.Unit or unit string (that can be converted to - a u.Unit object) represents some unit per solid angle. + Check if a given Unit or unit string (that can be converted to + a Unit) represents some unit per solid angle. If 'return_unit' + is True, then a Unit of the solid angle will be returned (or + None if no solid angle is present in the denominator). Parameters ---------- unit : str or u.Unit - Unit object or string representation of unit. + u.Unit object or string representation of unit. + return_unit : bool + If True, the u.Unit of the solid angle unit will + be returned (or None if unit is not a solid angle). Examples -------- @@ -139,7 +170,8 @@ def check_if_unit_is_per_solid_angle(unit): # this will take care of some formatting consistency like # turning something like Jy / (degree*degree) to Jy / deg**2 # and erg sr^1 to erg / sr - if isinstance(unit, u.core.Unit) or isinstance(unit, u.core.CompositeUnit): + if isinstance(unit, (u.core.Unit, u.core.CompositeUnit, + u.core.IrreducibleUnit)): unit_str = unit.to_string() elif isinstance(unit, str): unit = u.Unit(unit) @@ -148,16 +180,30 @@ def check_if_unit_is_per_solid_angle(unit): raise ValueError('Unit must be u.Unit, or string that can be converted into a u.Unit') if '/' in unit_str: - # might be comprised of several units in denom. + # input unit might be comprised of several units in denom. so check all. denom = unit_str.split('/')[-1].split() - # find all combos of one or two units, to catch cases where there are two different - # units of angle in the denom that might comprise a solid angle when multiplied. + # find all combos of one or two units, to catch cases where there are + # two different units of angle in the denom that might comprise a solid + # angle when multiplied. for i in [combo for length in (1, 2) for combo in itertools.combinations(denom, length)]: - # turn tuple of 1 or 2 units into a string, and turn that into a u.Unit to check type + # turn tuple of 1 or 2 units into a string, and turn that into a u.Unit + # to check type new_unit_str = ' '.join(i).translate(str.maketrans('', '', '()')) new_unit = u.Unit(new_unit_str) if new_unit.physical_type == 'solid angle': + if return_unit: # area units present and requested to be returned + return new_unit + return True # area units present but not requested to be returned + # square pixel should be considered a square angle unit + if new_unit == u.pix * u.pix: + if return_unit: + return new_unit return True + # in the case there are no area units, but return units were requested + if return_unit: + return None + + # and if there are no area units, and return units were NOT requested. return False diff --git a/jdaviz/tests/test_utils.py b/jdaviz/tests/test_utils.py index 86294bcd12..a099854b22 100644 --- a/jdaviz/tests/test_utils.py +++ b/jdaviz/tests/test_utils.py @@ -15,6 +15,8 @@ PHOTUTILS_LT_1_12_1 = not minversion(photutils, "1.12.1.dev") +PIX2 = u.pix * u.pix + def test_spec_sb_flux_conversion(): # Actual spectrum content does not matter, just the meta is used here. @@ -28,13 +30,21 @@ def test_spec_sb_flux_conversion(): assert_allclose(flux_conversion(values, u.Jy / u.sr, u.Jy, spec), [1, 2, 3]) assert_allclose(flux_conversion(values, u.Jy, u.Jy / u.sr, spec), [100, 200, 300]) + # conversions with eq pixels + assert_allclose(flux_conversion(values, u.Jy / PIX2, u.Jy, spec), [10, 20, 30]) + assert_allclose(flux_conversion(values, u.Jy, u.Jy / PIX2, spec), [10, 20, 30]) + # complex translation Jy / sr -> erg / (Angstrom s cm2 sr) targ = [2.99792458e-12, 1.49896229e-12, 9.99308193e-13] * (u.erg / (u.Angstrom * u.s * u.cm**2 * u.sr)) # noqa: E501 assert_allclose(flux_conversion(values, u.Jy / u.sr, u.erg / (u.Angstrom * u.s * u.cm**2 * u.sr), spec), targ.value) # noqa: E501 + # swap sr for pix2 and check conversion + assert_allclose(flux_conversion(values, u.Jy / PIX2, u.erg / (u.Angstrom * u.s * u.cm**2 * PIX2), spec), targ.value) # noqa: E501 # complex translation erg / (Angstrom s cm2 sr) -> Jy / sr targ = [3.33564095e+13, 2.66851276e+14, 9.00623057e+14] * (u.Jy / u.sr) assert_allclose(flux_conversion(values, u.erg / (u.Angstrom * u.s * u.cm**2 * u.sr), u.Jy / u.sr, spec), targ.value) # noqa: E501 + # swap sr for pix2 and check conversion + assert_allclose(flux_conversion(values, u.erg / (u.Angstrom * u.s * u.cm**2 * PIX2), u.Jy / PIX2, spec), targ.value) # noqa: E501 spectral_values = spec.spectral_axis spec_unit = u.MJy @@ -43,14 +53,17 @@ def test_spec_sb_flux_conversion(): # test spectrum when target unit in untranslatable unit list target_values = [5.03411657e-05, 2.01364663e-04, 4.53070491e-04] expected_units = (u.ph / (u.Hz * u.s * u.cm**2)) - returned_values, return_units, unit_flag = _indirect_conversion( - values=values, orig_units=(u.MJy), - targ_units=(u.ph / (u.s * u.cm**2 * u.Hz * u.sr)), # noqa - eqv=eqv, spec_unit=spec_unit, image_data=None - ) - assert_allclose(returned_values, target_values) - assert (return_units == expected_units) - assert (unit_flag == 'orig') + for solid_angle in [u.sr, u.pix*u.pix]: + returned_values, return_units, unit_flag = _indirect_conversion( + values=values, orig_units=(u.MJy), + targ_units=(u.ph / (u.s * u.cm**2 * u.Hz * solid_angle)), # noqa + eqv=eqv, + spec_unit=spec_unit, + image_data=None + ) + assert_allclose(returned_values, target_values) + assert (return_units == expected_units) + assert (unit_flag == 'orig') # test spectrum when original unit in untranslatable unit list target_values = [1., 2., 3.] diff --git a/jdaviz/utils.py b/jdaviz/utils.py index 389e274bf2..a1e7e93309 100644 --- a/jdaviz/utils.py +++ b/jdaviz/utils.py @@ -321,12 +321,20 @@ def standardize_roman_metadata(data_model): def indirect_units(): - return [ - u.erg / (u.s * u.cm**2 * u.Angstrom * u.sr), - u.erg / (u.s * u.cm**2 * u.Hz * u.sr), - u.ph / (u.Angstrom * u.s * u.cm**2 * u.sr), u.ph / (u.Angstrom * u.s * u.sr * u.cm**2), - u.ph / (u.s * u.cm**2 * u.Hz * u.sr) - ] + from jdaviz.core.validunits import supported_sq_angle_units + + units = [] + + for angle_unit in supported_sq_angle_units(): + units += [ + u.erg / (u.s * u.cm**2 * u.Angstrom * angle_unit), + u.erg / (u.s * u.cm**2 * u.Hz * angle_unit), + u.ph / (u.Angstrom * u.s * u.cm**2 * angle_unit), + u.ph / (u.Angstrom * u.s * angle_unit * u.cm**2), + u.ph / (u.s * u.cm**2 * u.Hz * angle_unit) + ] + + return units def flux_conversion(values, original_units, target_units, spec=None, eqv=None, slice=None): @@ -388,14 +396,15 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s eqv += u.spectral_density(slice) orig_units = u.Unit(original_units) - orig_bases = orig_units.bases targ_units = u.Unit(target_units) - targ_bases = targ_units.bases + + solid_angle_in_orig = check_if_unit_is_per_solid_angle(orig_units, return_unit=True) + 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 - (((u.sr in orig_bases) and (u.sr not in targ_bases)) or - ((u.sr not in orig_bases) and (u.sr in targ_bases)))): + (((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. # App-wide original data units are used for conversion, original and # target_units dictate the conversion to take place. @@ -413,6 +422,12 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s eqv_in = fac eqv += _eqv_pixar_sr(np.array(eqv_in)) + # may need equivalencies between flux and flux per square pixel + eqv += _eqv_flux_to_sb_pixel() + + # when angle<>pixel translations are enabled + # eqv += _eqv_sb_per_pixel_to_per_angle(u.Jy) + # indirect units cannot be directly converted, and require # additional conversions to reach the desired end unit. # if spec_unit in [original_units, target_units]: @@ -436,41 +451,62 @@ def flux_conversion(values, original_units, target_units, spec=None, eqv=None, s values=values, orig_units=orig_units, targ_units=targ_units, eqv=eqv, image_data=image_data ) + elif solid_angle_in_orig == solid_angle_in_targ == u.pix * u.pix: + # in the case where we have 2 SBs per solid pixel that need + # u.spectral_density equivalency, they can't be directly converted + # for whatever reason (i.e 'Jy / pix2' and 'erg / (Angstrom s cm2 pix2)' + # are not convertible). In this case, multiply out the factor of pix2 for + # conversion (same kind of thing _indirect_conversion is + # doing but we already know the exact angle units. + orig_units *= u.pix * u.pix + targ_units *= u.pix * u.pix return (values * orig_units).to_value(targ_units, equivalencies=eqv) def _indirect_conversion(values, orig_units, targ_units, eqv, spec_unit=None, image_data=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 + # assume pix2 now if there are none? or use the display units? + + solid_angle_in_orig = check_if_unit_is_per_solid_angle(orig_units, return_unit=True) + solid_angle_in_targ = check_if_unit_is_per_solid_angle(targ_units, return_unit=True) + if spec_unit is not None: + solid_angle_in_spec = check_if_unit_is_per_solid_angle(spec_unit, return_unit=True) + else: + solid_angle_in_spec = None + # indirect units cannot be directly converted, and require # additional conversions to reach the desired end unit. - if (spec_unit and spec_unit in [orig_units, targ_units] - and not check_if_unit_is_per_solid_angle(spec_unit)): + if (spec_unit and spec_unit in [orig_units, targ_units] and not solid_angle_in_spec): if u.Unit(targ_units) in indirect_units(): - temp_targ = targ_units * u.sr + temp_targ = targ_units * solid_angle_in_targ values = (values * orig_units).to_value(temp_targ, equivalencies=eqv) orig_units = u.Unit(temp_targ) return values, orig_units, 'orig' elif u.Unit(orig_units) in indirect_units(): - temp_orig = orig_units * u.sr + temp_orig = orig_units * solid_angle_in_orig values = (values * orig_units).to_value(temp_orig, equivalencies=eqv) targ_units = u.Unit(temp_orig) return values, targ_units, 'targ' return values, targ_units, 'targ' - elif image_data or (spec_unit and check_if_unit_is_per_solid_angle(spec_unit)): - if not check_if_unit_is_per_solid_angle(targ_units): - targ_units /= u.sr + elif image_data 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 (u.Unit(orig_units) in indirect_units())): # SB -> Flux -> Flux -> SB - temp_orig = orig_units * u.sr - temp_targ = targ_units * u.sr + temp_orig = orig_units * solid_angle_in_orig + temp_targ = targ_units * solid_angle_in_targ # Convert Surface Brightness to Flux, then Flux to Flux values = (values * orig_units).to_value(temp_orig, equivalencies=eqv) values = (values * temp_orig).to_value(temp_targ, equivalencies=eqv) + # Lastly a Flux to Surface Brightness translation in the return statement orig_units = temp_targ @@ -480,6 +516,11 @@ def _indirect_conversion(values, orig_units, targ_units, eqv, def _eqv_pixar_sr(pixar_sr): + """ + Return Equivalencies to convert from flux to flux per solid + angle (aka surface brightness) using scale ratio `pixar_sr` + (steradians per pixel). + """ def converter_flux(x): # Surface Brightness -> Flux return x * pixar_sr @@ -490,10 +531,60 @@ def iconverter_flux(x): # Flux -> Surface Brightness (u.MJy / u.sr, u.MJy, converter_flux, iconverter_flux), (u.erg / (u.s * u.cm**2 * u.Angstrom * u.sr), u.erg / (u.s * u.cm**2 * u.Angstrom), converter_flux, iconverter_flux), # noqa (u.ph / (u.Angstrom * u.s * u.cm**2 * u.sr), u.ph / (u.Angstrom * u.s * u.cm**2), converter_flux, iconverter_flux), # noqa - (u.ph / (u.Hz * u.s * u.cm**2 * u.sr), u.ph / (u.Hz * u.s * u.cm**2), converter_flux, iconverter_flux) # noqa + (u.ph / (u.Hz * u.s * u.cm**2 * u.sr), u.ph / (u.Hz * u.s * u.cm**2), converter_flux, iconverter_flux), # noqa + (u.ct / u.sr, u.ct, converter_flux, iconverter_flux) # noqa ] +def _eqv_flux_to_sb_pixel(): + """ + Returns an Equivalency between `flux_unit` and `flux_unit`/pix**2. This + allows conversion between flux and flux-per-square-pixel surface brightness + e.g MJy <> MJy / pix2 + """ + + pix2 = u.pix * u.pix + + # generate an equivalency for each flux type that would need + # another equivalency for converting to/from + flux_units = [u.MJy, u.erg / (u.s * u.cm**2 * u.Angstrom), + u.ph / (u.Angstrom * u.s * u.cm**2), + u.ph / (u.Hz * u.s * u.cm**2)] + equiv = [] + for flux_unit in flux_units: + equiv.append((flux_unit, flux_unit / pix2, lambda x: x, lambda x: x)) + + return equiv + + +def _eqv_sb_per_pixel_to_per_angle(flux_unit, scale_factor=1): + """ + Returns an equivalency between `flux_unit` per square pixel and + `flux_unit` per solid angle to be able to compare and convert between units + like Jy/pix**2 and Jy/sr. The scale factor is assumed to be in steradians, + to follow the convention of the PIXAR_SR keyword. + Note: + To allow conversions between units like 'ph / (Hz s cm2 sr)' and + MJy / pix2, which would require this equivalency as well as u.spectral_density, + these CAN'T be combined when converting like: + equivalencies=u.spectral_density(1 * u.m) + _eqv_sb_per_pixel_to_per_angle(u.Jy) + So additional logic is needed to compare units that need both equivalencies + (one solution being creating this equivalency for each equivalent flux-type.) + + """ + pix2 = u.pix * u.pix + + # the two types of units we want to define a conversion between + flux_solid_ang = flux_unit / u.sr + flux_sq_pix = flux_unit / pix2 + + pix_to_solid_angle_equiv = [(flux_solid_ang, flux_sq_pix, + lambda x: x * scale_factor, + lambda x: x / scale_factor)] + + return pix_to_solid_angle_equiv + + def spectral_axis_conversion(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)