From ca3e40e71f3b2fdfca2dea7eb87bc3a15990f54f Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Fri, 6 Sep 2024 01:51:26 -0400 Subject: [PATCH] checkpoint --- jdaviz/app.py | 11 +- .../plugins/moment_maps/moment_maps.py | 60 +++------ .../moment_maps/tests/test_moment_maps.py | 43 ++++-- jdaviz/configs/cubeviz/plugins/parsers.py | 123 ++++++++++++++++-- .../spectral_extraction.py | 15 ++- .../tests/test_spectral_extraction.py | 5 + .../plugins/tests/test_cubeviz_aperphot.py | 83 +++++++++--- .../cubeviz/plugins/tests/test_parsers.py | 8 +- .../cubeviz/plugins/tests/test_tools.py | 22 +++- .../markers/tests/test_markers_plugin.py | 22 ++-- .../model_fitting/tests/test_fitting.py | 27 +++- jdaviz/configs/default/plugins/viewers.py | 56 +++++--- .../aper_phot_simple/aper_phot_simple.py | 93 +++++++++---- .../aper_phot_simple/aper_phot_simple.vue | 3 +- .../unit_conversion/unit_conversion.py | 38 ++++-- jdaviz/conftest.py | 5 + jdaviz/core/marks.py | 5 +- jdaviz/core/validunits.py | 109 +++++++++++----- jdaviz/utils.py | 119 ++++++++++++++--- 19 files changed, 630 insertions(+), 217 deletions(-) diff --git a/jdaviz/app.py b/jdaviz/app.py index 5ba4f47752..e022d92bb4 100644 --- a/jdaviz/app.py +++ b/jdaviz/app.py @@ -86,7 +86,16 @@ def equivalent_units(self, data, cid, units): 'erg / (s sr cm2)', 'erg / (Hz s sr cm2)', 'erg / (Angstrom s sr cm2)', 'ph / (Angstrom s sr cm2)', 'ph / (Hz s sr cm2)' - ]) + ] + + [ + '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)' + ] + + ) else: # spectral axis # prefer Hz over Bq and um over micron exclude = {'Bq', 'micron'} diff --git a/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py b/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py index 2c71128fd9..706f3c0142 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_flux_or_sb = self.output_unit_items[0]['label'] - - unit_dict = {orig_flux_or_sb: "", + 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] @@ -197,36 +193,11 @@ def _set_data_units(self, event={}): 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_flux_or_sb] = 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_flux_or_sb] = 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_flux_or_sb] = 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_flux_or_sb]) - 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,21 +314,20 @@ 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 - # convert unit string to u.Unit so moment map data can be converted - flux_or_sb_display_unit = u.Unit(flux_sb_unit) + # get display units for moment 0 based on unit conversion plugin selection + moment_0_display_unit = self.output_unit_items[0]['unit_str'] + + # convert unit string to Unit so moment map data can be converted + flux_or_sb_display_unit = u.Unit(moment_0_display_unit) + + # account for extra wavelength factor, depending on specutils version if SPECUTILS_LT_1_15_1: moment_new_unit = flux_or_sb_display_unit else: moment_new_unit = flux_or_sb_display_unit * self.spectrum_viewer.state.x_display_unit # noqa: E501 + + print(self.moment.unit, moment_new_unit) 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 859ec48b37..0a15a0de2e 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 @@ -14,10 +14,21 @@ from jdaviz.configs.cubeviz.plugins.moment_maps.moment_maps import SPECUTILS_LT_1_15_1 -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 @@ -45,12 +56,15 @@ 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): + if SPECUTILS_LT_1_15_1: moment_unit = "Jy" moment_value_str = "+8.00000e+00" @@ -58,10 +72,21 @@ def test_moment_calculation(cubeviz_helper, spectrum1d_cube, 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 = moment_unit # 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. @@ -73,7 +98,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() @@ -93,13 +120,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)") diff --git a/jdaviz/configs/cubeviz/plugins/parsers.py b/jdaviz/configs/cubeviz/plugins/parsers.py index d9b2484a11..2deac9eae7 100644 --- a/jdaviz/configs/cubeviz/plugins/parsers.py +++ b/jdaviz/configs/cubeviz/plugins/parsers.py @@ -8,9 +8,11 @@ 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.configs.imviz.plugins.parsers import prep_data_layer_as_dq 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 +179,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 +295,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 +311,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,11 +358,7 @@ 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 - data = app.data_collection[-1] prep_data_layer_as_dq(data) @@ -448,6 +461,12 @@ 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": + if 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 +493,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 +520,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 +552,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 +574,71 @@ 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 `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 adjusting the units of spectrum.flux + and (if present) spectrum.uncertainty, and creating a new `Spectrum1D` + object with the modified flux and uncertainty. + + Parameters + ---------- + spectrum : Spectrum1D + A `Spectrum1D` object containing flux data, which is assumed to be in + flux units without any angular component in the denominator. + + Returns + ------- + Spectrum1D + A new `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') + old_flux_unit = flux.unit + unitless_flux = flux.value + flux = unitless_flux * old_flux_unit / (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. + old_uncerts_unit = old_uncerts.unit + unitless_uncerts = uncerts.array + uncerts = unitless_uncerts * old_uncerts_unit / (u.pix * u.pix) + uncerts = StdDevUncertainty(uncerts) + + # 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 signiture 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 b82acec4f7..08ccdebe97 100644 --- a/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py +++ b/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py @@ -24,6 +24,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 @@ -503,9 +504,17 @@ 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) + print('sq angle unit', sq_angle_unit) + 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 295432051e..7ef1a4a125 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 c600ee01c9..622dabe7c5 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,14 +124,20 @@ 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) -def test_cubeviz_aperphot_cube_orig_flux_mjysr(cubeviz_helper, spectrum1d_cube_custom_fluxunit): - cube = spectrum1d_cube_custom_fluxunit(fluxunit=u.MJy / u.sr) +@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): + 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) @@ -132,9 +149,29 @@ def test_cubeviz_aperphot_cube_orig_flux_mjysr(cubeviz_helper, spectrum1d_cube_c plg.aperture_selected = "Subset 1" plg.background_selected = "Subset 2" - # Make sure per steradian is handled properly. - assert_allclose(plg.pixel_area, 0.01) - assert_allclose(plg.flux_scaling, 0.003631) + # 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 + assert not plg.disable_pixarea_input # for /sr cubes, this should be editable + + 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 + assert plg.disable_pixarea_input plg.vue_do_aper_phot() row = cubeviz_helper.get_aperture_photometry_results()[0] @@ -142,11 +179,19 @@ def test_cubeviz_aperphot_cube_orig_flux_mjysr(cubeviz_helper, spectrum1d_cube_c # Basically, we should recover the input rectangle here, minus background. assert_allclose(row["xcenter"], 3 * u.pix) assert_allclose(row["ycenter"], 1 * u.pix) - assert_allclose(row["sum"], 1.1752215e-12 * u.MJy) # (15 - 10) MJy/sr x 2.3504431e-13 sr + # (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)) - assert_allclose(row["pixarea_tot"], 2.350443053909789e-13 * u.sr) - assert_allclose(row["aperture_sum_mag"], 23.72476627732448 * u.mag) - assert_allclose(row["mean"], 5 * (u.MJy / u.sr)) + + # 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) diff --git a/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py b/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py index 8946cc283e..557bf36882 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,7 +152,7 @@ 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' def test_spectrum1d_parse(spectrum1d, cubeviz_helper): 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..eeb97618fa 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}', '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}', '', '') 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/tests/test_fitting.py b/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py index cbddd12c59..00eefdf90d 100644 --- a/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py +++ b/jdaviz/configs/default/plugins/model_fitting/tests/test_fitting.py @@ -76,12 +76,22 @@ def test_model_ids(cubeviz_helper, spectral_cube_wcs): 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.*') @@ -90,14 +100,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)) @@ -223,6 +234,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) @@ -273,7 +286,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 == 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], diff --git a/jdaviz/configs/default/plugins/viewers.py b/jdaviz/configs/default/plugins/viewers.py index 5bf05593e6..4957d34249 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,26 +756,40 @@ 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: - # default to Flux Density for flux density or uncaught types - flux_unit_type = "Flux density" + 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 + # defaulty to the catchall 'flux density' label + flux_unit_type = None + + if solid_angle_unit is not None: + locally_defined_sb_units = [un / solid_angle_unit for un in locally_defined_flux_units] + + # create an equivalency for each flux unit for flux <> flux/pix2. + # for similar reasons to the 'untranslatable units' issue, custom equivalencies can't be combined, + # so a workaround is creating an equivalency for each flux that may need an additional equiv. + angle_to_pixel_equivs = [_eqv_sb_per_pixel_to_per_angle(unit) for unit in locally_defined_flux_units] + + if any([y_unit.is_equivalent(unit, angle_to_pixel_equivs[i]) for i, unit in enumerate(locally_defined_sb_units)]): # noqa + flux_unit_type = "Surface Brightness" + + if flux_unit_type is None: # if a label of SB or Flux density wasn't determined + if any(y_unit.is_equivalent(unit) for unit in locally_defined_flux_units): + flux_unit_type = 'Flux' + elif 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" + else: + # default to Flux Density for flux density or uncaught types + flux_unit_type = "Flux density" # Set x axes labels for the spectrum viewer x_disp_unit = self.state.x_display_unit 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 7a9768a6a9..251f074877 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, @@ -69,7 +71,9 @@ class SimpleAperturePhotometry(PluginTemplateMixin, ApertureSubsetSelectMixin, cube_slice = Unicode("").tag(sync=True) is_cube = Bool(False).tag(sync=True) display_flux_or_sb_unit = Unicode("").tag(sync=True) + display_solid_angle_unit = Unicode("").tag(sync=True) flux_scaling_display_unit = Unicode("").tag(sync=True) + disable_pixarea_input = Bool(False).tag(sync=True) # if flux/pix2 def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) @@ -108,16 +112,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, @@ -191,6 +202,7 @@ def _on_display_units_changed(self, event={}): new_unit = u.Unit(self.display_flux_or_sb_unit) bg = self.background_value * prev_unit + self.background_value = bg.to_value( new_unit, u.spectral_density(self._cube_wave)) @@ -212,11 +224,16 @@ def _set_display_unit_of_selected_dataset(self): surface brightness. """ + # the refactoring done in 3144 will have the response to + # GlobalDisplayUnitChanged handled correctly. for now, use + # _get_display_unit + if not self.dataset_selected or not self.aperture_selected: self.display_flux_or_sb_unit = '' self.flux_scaling_display_unit = '' return +<<<<<<< HEAD data = self.dataset.selected_dc_item if isinstance(data, list): data = data[0] @@ -228,20 +245,32 @@ def _set_display_unit_of_selected_dataset(self): flux_or_sb = 'sb' else: flux_or_sb = 'flux' +======= + disp_unit = self.app._get_display_unit('sb') # all cubes are in sb + self.display_flux_or_sb_unit = disp_unit +>>>>>>> efb22333 (added ap phot changes, all cubeviz plugin test passing) - disp_unit = self.app._get_display_unit(flux_or_sb) + # get angle componant of surface brightness + # note: could add 'axis=angle' when cleaning this code up to avoid repeating this - self.display_flux_or_sb_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_flux_or_sb_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.disable_pixarea_input = True + self.pixel_area = 1.0 def _get_defaults_from_metadata(self, dataset=None): defaults = {} @@ -349,7 +378,7 @@ def _dataset_selected_changed(self, event={}): # get correct display unit for newly selected dataset if self.config == 'cubeviz': - # set display_flux_or_sb_unit and flux_scaling_display_unit + # sets display_flux_or_sb_unit and flux_scaling_display_unit traitlets self._set_display_unit_of_selected_dataset() # auto-populate background, if applicable. @@ -501,7 +530,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 @@ -584,8 +613,9 @@ def calculate_photometry(self, dataset=None, aperture=None, background=None, background_value = (background_value * display_unit).to_value( img_unit, u.spectral_density(self._cube_wave)) else: - bg_reg = self.aperture._get_spatial_region(subset=background if background is not None else self.background.selected, # noqa - dataset=dataset if dataset is not None else self.dataset.selected) # noqa + sub = background if background is not None else self.background.selected + dat = dataset if dataset is not None else self.dataset.selected + bg_reg = self.aperture._get_spatial_region(subset=sub, dataset=dat) background_value = self._calc_background_median(bg_reg, data=data) # cubeviz: computed background median will be in display units, @@ -620,7 +650,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: @@ -663,9 +694,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) @@ -689,11 +722,18 @@ 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. - # 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)) + # convert pixarea, which is in arcsec2/pix2 to the display solid angle unit / pix2 + 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) + phot_table['sum'] = [rawsum * pixarea_fac] else: pixarea_fac = None @@ -1223,7 +1263,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 7391f01d86..f4432d678e 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 /> - + + Flux - if check_if_unit_is_per_solid_angle(spec_units) and flux_or_sb == 'Flux': - spec_units *= u.sr + if spec_axis_ang_unit and flux_or_sb == 'Flux': + print(f'translating SB to flux. spec_units = {spec_units}, angle={selected_display_solid_angle_unit}') + spec_units *= selected_display_solid_angle_unit + print('new spec_units = ', spec_units) # 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 flux_or_sb == 'Surface Brightness'): - spec_units /= u.sr + elif (not spec_axis_ang_unit and flux_or_sb == 'Surface Brightness'): + print(f'translating flux to sb. spec_units = {spec_units}, angle={selected_display_solid_angle_unit}') + spec_units /= selected_display_solid_angle_unit + print('new spec_units = ', spec_units) # update display units self.spectrum_viewer.state.y_display_unit = str(spec_units) + # entered the translator when we shouldn't translate else: return @@ -351,7 +369,7 @@ def _translate(self, flux_or_sb=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: diff --git a/jdaviz/conftest.py b/jdaviz/conftest.py index 0be26cb176..89b42e7db8 100644 --- a/jdaviz/conftest.py +++ b/jdaviz/conftest.py @@ -305,6 +305,11 @@ def spectrum1d_cube_fluxunit_jy_per_steradian(): return _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy/u.sr, shape=(10, 4, 5), 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 af21d275b2..0ff1c56d80 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, @@ -116,6 +116,9 @@ def set_y_unit(self, unit=None): 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) + + # for flux <> flux/pix2 + eqv += _eqv_flux_to_sb_pixel(u.Jy) # extend to untranslatable, Jy is placeholder 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..b20b5825a3 100644 --- a/jdaviz/core/validunits.py +++ b/jdaviz/core/validunits.py @@ -1,9 +1,11 @@ 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 +75,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,44 +95,70 @@ 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. + + """ -def check_if_unit_is_per_solid_angle(unit): + 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] + + 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. + a u.Unit object) represents some unit per solid angle. If 'return_unit' + is True, then a u.Unit of the solid angle in `unit` 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. + 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 -------- >>> check_if_unit_is_per_solid_angle('erg / (s cm^2 sr)') True + >>> check_if_unit_is_per_solid_angle('erg / (s cm^2 sr)', return_unit=True) + Unit("sr") >>> check_if_unit_is_per_solid_angle('erg / s cm^2') False >>> check_if_unit_is_per_solid_angle('Jy * sr^-1') @@ -141,6 +172,8 @@ def check_if_unit_is_per_solid_angle(unit): # 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, u.core.IrreducibleUnit): # u.count + unit_str = unit.to_string() elif isinstance(unit, str): unit = u.Unit(unit) unit_str = unit.to_string() @@ -148,16 +181,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/utils.py b/jdaviz/utils.py index 389e274bf2..ddd84f2265 100644 --- a/jdaviz/utils.py +++ b/jdaviz/utils.py @@ -320,13 +320,28 @@ 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) +# ] + 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 +403,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) + solid_angle_in_targ = check_if_unit_is_per_solid_angle(targ_units) # 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)))): + if (('_pixel_scale_factor' in spec.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. # App-wide original data units are used for conversion, original and # target_units dictate the conversion to take place. @@ -413,6 +429,14 @@ 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 + # NOTE need to extend this to other flux types that aren't equiv with u.Jy, + # this is a placeholder + eqv += _eqv_flux_to_sb_pixel(u.Jy) + + # 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]: @@ -442,35 +466,48 @@ 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): + + print(f'in _indirect_conversion, orig_units={orig_units}, targ_units={targ_units}') + + # make these an input parameter since they're already determined + # here for now until i make sure this function isn't called elsewhere + 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) + solid_angle_in_spec = check_if_unit_is_per_solid_angle(spec_unit, return_unit=True) + # 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 u.Unit(targ_units) in indirect_units(): - temp_targ = targ_units * u.sr + temp_targ = targ_units * solid_angle_in_targ + print('temp_targ', temp_targ) + print('solid angle in targ', solid_angle_in_targ) + print('orig_units', orig_units) 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 @@ -478,8 +515,15 @@ def _indirect_conversion(values, orig_units, targ_units, eqv, return values, orig_units +# Note: should unify how these custom equivs are written, either they should take in any flux unit +# or return multiple equivalencies for all non-equivalent flux unit types. 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 @@ -493,6 +537,45 @@ def iconverter_flux(x): # Flux -> Surface Brightness (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 ] +def _eqv_flux_to_sb_pixel(flux_unit): + """ + 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 + equiv = [(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 = u.Equivalency([(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)