Skip to content
New issue

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

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

Already on GitHub? # to your account

Correct sum units for per steradian #2873

Merged
merged 12 commits into from
Jun 26, 2024

Conversation

javerbukh
Copy link
Contributor

Description

This pull request sets the extracted spectrum to the correct units when using the sum function with flux units in per steradian.

Change log entry

  • Is a change log needed? If yes, is it added to CHANGES.rst? If you want to avoid merge conflicts,
    list the proposed change log here for review and add to CHANGES.rst before merge. If no, maintainer
    should add a no-changelog-entry-needed label.
  • Set correct sum units when flux units are in per steradian.

Checklist for package maintainer(s)

This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.

  • Are two approvals required? Branch protection rule does not check for the second approval. If a second approval is not necessary, please apply the trivial label.
  • Do the proposed changes actually accomplish desired goals? Also manually run the affected example notebooks, if necessary.
  • Do the proposed changes follow the STScI Style Guides?
  • Are tests added/updated as required? If so, do they follow the STScI Style Guides?
  • Are docs added/updated as required? If so, do they follow the STScI Style Guides?
  • Did the CI pass? If not, are the failures related?
  • Is a milestone set? Set this to bugfix milestone if this is a bug fix and needs to be released ASAP; otherwise, set this to the next major release milestone. Bugfix milestone also needs an accompanying backport label.
  • After merge, any internal documentations need updating (e.g., JIRA, Innerspace)?

@github-actions github-actions bot added cubeviz testing plugin Label for plugins common to multiple configurations labels May 13, 2024
@javerbukh javerbukh added this to the 3.11 milestone May 13, 2024
Comment on lines 378 to 439
new_unit = collapsed_nddata.unit * astropy.units.sr
collapsed_nddata.unit = new_unit
if hasattr(collapsed_nddata.uncertainty, 'unit'):
new_uncert_unit = collapsed_nddata.uncertainty.unit * astropy.units.sr
collapsed_nddata.uncertainty.unit = new_uncert_unit
Copy link
Contributor

Choose a reason for hiding this comment

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

This conversion needs a multiplication with an area term that has units [sr/pixel]. The pixels have units like [Jy / sr]. If we want to sum over an area and remove the area unit, we need to multiply (Jy / sr) * (sr / pixel). The second term is stored in the pixar_sr keyword, and you can get it with:

self.spectral_cube.meta.get('PIXAR_SR', 1.0)

Copy link
Contributor

Choose a reason for hiding this comment

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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Does the latest commit look correct?

Copy link
Member

Choose a reason for hiding this comment

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

you need the area of the aperture itself. You can just refer to self.aperture_area_along_spectral directly to get the area (in pixels) and multiply by the PIXAR_SR to get it in sr.

Copy link
Contributor Author

@javerbukh javerbukh May 29, 2024

Choose a reason for hiding this comment

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

self.aperture_area_along_spectral is an array, so is the area one of the elements of the array or the elements multiplied together?

Edit: I hadn't had enough coffee, the test data had a spectral axis with a length of 2 and that tripped me up 🤦‍♂️ I understand now.

Copy link
Member

Choose a reason for hiding this comment

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

its an area for each slice (since the aperture could be a cone)

Copy link
Contributor Author

@javerbukh javerbukh May 29, 2024

Choose a reason for hiding this comment

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

I see that the wavelength_dependent parameter is not used in _extract_from_aperture, does that mean main no longer has that functionality?

Copy link
Member

Choose a reason for hiding this comment

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

It can probably be removed from the method signature, yes. Wavelength dependence is now being handled by passing reference_spectral_value to get_mask in the aperture_weight_mask and bg_weight_mask properties, respectively (although please feel free to test to make sure its working as expected).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Hopefully I'm on the right track now. I added a check of the collapsed flux to the test so you can see if the values look correct. The flux from the original collapsed nddata sum array have been multiplied by the self.aperture_area_along_spectral values (in the test case they are 20).

@javerbukh javerbukh force-pushed the spec-extract-sum-units branch from 71b83ef to 8e0d989 Compare May 29, 2024 13:48
axis=spatial_axes, **kwargs
) # returns an NDDataArray
# Remove per steradian denominator
if astropy.units.sr in collapsed_nddata.unit.bases:
Copy link
Contributor

Choose a reason for hiding this comment

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

If you want to check that the unit is per solid angle (rather than just checking if the unit contains steradian), you can use core.validunits.check_if_unit_is_per_solid_angle. that would require some additional logic, which you could add to that function, to return what the solid angle unit is to multiply by.

Comment on lines 437 to 440
collapsed_nddata.unit *= aperture_area.unit
collapsed_nddata.data[:] *= aperture_area
# Update uncertainty units as well
collapsed_nddata.uncertainty[:] *= aperture_area
Copy link
Contributor

Choose a reason for hiding this comment

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

Let's use NDData.multiply() instead (docs).

@javerbukh
Copy link
Contributor Author

Test failures are present on main but made visible in this PR. The core issue seems to be the _handle_display_units() does not handle the surface brightness to flux conversion, and all plugins that call _get_data() eventually will hit this issue. To trigger, follow these steps:

  1. Run CubevizExample notebook to data load step

2.Run the following lines of code

uc = cubeviz.plugins['Unit Conversion']._obj
uc.show_translator = True
uc.flux_or_sb = 'Flux'
la = cubeviz.plugins['Line Analysis']
la.get_results()
  1. The traceback should look like this:
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
Cell In[5], line 2
      1 la = cubeviz.plugins['Line Analysis']
----> 2 la.get_results()

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py:219](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py#line=218), in LineAnalysis.get_results(self)
    215     valid, spec_range, subset_range = self._check_dataset_spectral_subset_valid(return_ranges=True)  # noqa
    216     raise ValueError(f"spectral subset '{self.spectral_subset.selected}' {subset_range}"
    217                      f" is outside data range of '{self.dataset.selected}' {spec_range}")
--> 219 self._calculate_statistics()
    220 return self.results

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py:318](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py#line=317), in skip_if_no_updates_since_last_active.<locals>.decorator.<locals>.wrapper(self, msg)
    316 if meth.__name__ not in self._methods_skip_since_last_active:
    317     self._methods_skip_since_last_active.append(meth.__name__)
--> 318 ret_ = meth(self, msg)
    319 if ret_ is False:
    320     self._methods_skip_since_last_active.remove(meth.__name__)

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py:349](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py#line=348), in with_spinner.<locals>.decorator.<locals>.wrapper(self, *args, **kwargs)
    347 setattr(self, spinner_traitlet, True)
    348 try:
--> 349     ret_ = meth(self, *args, **kwargs)
    350 except Exception:
    351     setattr(self, spinner_traitlet, False)

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py:263](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py#line=262), in LineAnalysis._calculate_statistics(self, msg)
    260     self.update_results(None)
    261     return
--> 263 spectrum, continuum, spec_subtracted = self._get_continuum(self.dataset,
    264                                                            self.spectral_subset,
    265                                                            update_marks=True)
    266 if spectrum is None:
    267     self.update_results(None)

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py:2950](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py#line=2949), in SpectralContinuumMixin._get_continuum(self, dataset, spectral_subset, update_marks, per_pixel)
   2947     full_spectrum = self.app._jdaviz_helper.get_data(self.dataset.selected,
   2948                                                      use_display_units=True)
   2949 else:
-> 2950     full_spectrum = dataset.get_selected_spectrum(use_display_units=True)
   2952 if full_spectrum is None or self.continuum_width == "":
   2953     self._update_continuum_marks()

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py:3437](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/template_mixin.py#line=3436), in DatasetSelect.get_selected_spectrum(self, use_display_units)
   3433     sp = spec_extract._extract_in_new_instance(self.selected,
   3434                                                function=self._spectral_extraction_function,
   3435                                                add_data=False)
   3436     return self.plugin._specviz_helper._handle_display_units(sp, use_display_units)
-> 3437 return self.plugin._specviz_helper.get_data(data_label=self.selected,
   3438                                             use_display_units=use_display_units)

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/helper.py:313](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/configs/specviz/helper.py#line=312), in Specviz.get_data(self, data_label, spectral_subset, cls, use_display_units)
    290 def get_data(self, data_label=None, spectral_subset=None, cls=None,
    291              use_display_units=False):
    292     """
    293     Returns data with name equal to data_label of type cls with subsets applied from
    294     spectral_subset.
   (...)
    311 
    312     """
--> 313     return self._get_data(data_label=data_label, spectral_subset=spectral_subset,
    314                           cls=cls, use_display_units=use_display_units)

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/helpers.py:548](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/helpers.py#line=547), in ConfigHelper._get_data(self, data_label, spatial_subset, spectral_subset, mask_subset, cls, use_display_units)
    545     else:
    546         data = data.get_object(cls=cls, **object_kwargs)
--> 548     return self._handle_display_units(data, use_display_units)
    550 if not cls and spatial_subset:
    551     raise AttributeError(f"A valid cls must be provided to"
    552                          f" apply subset {spatial_subset} to data. "
    553                          f"Instead, {cls} was given.")

File [~/Documents/jdaviz_dev/jdaviz/jdaviz/core/helpers.py:477](http://localhost:8889/~/Documents/jdaviz_dev/jdaviz/jdaviz/core/helpers.py#line=476), in ConfigHelper._handle_display_units(self, data, use_display_units)
    473     uncertainty.unit = data.flux.unit
    474 if hasattr(uncertainty, 'represent_as'):
    475     new_uncert = uncertainty.represent_as(
    476         StdDevUncertainty
--> 477     ).quantity.to(flux_unit)
    478 else:
    479     # if not specified as NDUncertainty, assume stddev:
    480     new_uncert = uncertainty.quantity.to(flux_unit)

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/quantity.py:938](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/quantity.py#line=937), in Quantity.to(self, unit, equivalencies, copy)
    934 unit = Unit(unit)
    935 if copy:
    936     # Avoid using to_value to ensure that we make a copy. We also
    937     # don't want to slow down this method (esp. the scalar case).
--> 938     value = self._to_value(unit, equivalencies)
    939 else:
    940     # to_value only copies if necessary
    941     value = self.to_value(unit, equivalencies)

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/quantity.py:891](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/quantity.py#line=890), in Quantity._to_value(self, unit, equivalencies)
    888     equivalencies = self._equivalencies
    889 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    890     # Standard path, let unit to do work.
--> 891     return self.unit.to(
    892         unit, self.view(np.ndarray), equivalencies=equivalencies
    893     )
    895 else:
    896     # The .to() method of a simple unit cannot convert a structured
    897     # dtype, so we work around it, by recursing.
    898     # TODO: deprecate this?
    899     # Convert simple to Structured on initialization?
    900     result = np.empty_like(self.view(np.ndarray))

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py:1195](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py#line=1194), in UnitBase.to(self, other, value, equivalencies)
   1193     return UNITY
   1194 else:
-> 1195     return self._get_converter(Unit(other), equivalencies)(value)

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py:1124](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py#line=1123), in UnitBase._get_converter(self, other, equivalencies)
   1121             else:
   1122                 return lambda v: b(converter(v))
-> 1124 raise exc

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py:1107](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py#line=1106), in UnitBase._get_converter(self, other, equivalencies)
   1105 # if that doesn't work, maybe we can do it with equivalencies?
   1106 try:
-> 1107     return self._apply_equivalencies(
   1108         self, other, self._normalize_equivalencies(equivalencies)
   1109     )
   1110 except UnitsError as exc:
   1111     # Last hope: maybe other knows how to do it?
   1112     # We assume the equivalencies have the unit itself as first item.
   1113     # TODO: maybe better for other to have a `_back_converter` method?
   1114     if hasattr(other, "equivalencies"):

File [~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py:1085](http://localhost:8889/~/miniconda3/envs/jdaviz-dev-3.9/lib/python3.10/site-packages/astropy/units/core.py#line=1084), in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
   1082 unit_str = get_err_str(unit)
   1083 other_str = get_err_str(other)
-> 1085 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: 'MJy [/](http://localhost:8889/) sr' and 'MJy' (spectral flux density) are not convertible

Copy link
Contributor

@bmorris3 bmorris3 left a comment

Choose a reason for hiding this comment

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

This looks great.

if astropy.units.sr in collapsed_nddata.unit.bases:
aperture_area = (self.aperture_area_along_spectral
* self.spectral_cube.meta.get('PIXAR_SR', 1.0) * u.sr)
collapsed_nddata = collapsed_nddata.multiply(aperture_area,
Copy link
Contributor

Choose a reason for hiding this comment

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

Lovely 👏🏻


def test_spectral_extraction_with_correct_sum_units(cubeviz_helper,
spectrum1d_cube_fluxunit_jy_per_steradian):
cubeviz_helper.load_data(spectrum1d_cube_fluxunit_jy_per_steradian)
Copy link
Contributor

Choose a reason for hiding this comment

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

This is not strictly necessary, but good for future quality of life. Could we add an assert after this line to ensure that the loaded cube's flux unit is Jy/sr? That would make it really clear what this test is checking for later on.

@@ -214,7 +214,7 @@ def _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy, shape=(2, 2, 4), with_u
"TELESCOP": "JWST", "BUNIT": fluxunit.to_string(), "PIXAR_A2": 0.01}
w = WCS(wcs_dict)
if with_uncerts:
uncert = StdDevUncertainty(np.abs(np.random.normal(flux) * u.Jy))
uncert = StdDevUncertainty(np.abs(np.random.normal(flux) * fluxunit))
Copy link
Contributor

Choose a reason for hiding this comment

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

Good catch. 😅

@javerbukh javerbukh force-pushed the spec-extract-sum-units branch from 278c7d8 to b6eb788 Compare June 20, 2024 11:46
@javerbukh javerbukh force-pushed the spec-extract-sum-units branch from b6eb788 to d6e370a Compare June 20, 2024 11:57
Copy link

codecov bot commented Jun 20, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 88.79%. Comparing base (55e6f6a) to head (87879f0).
Report is 148 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff           @@
##             main    #2873   +/-   ##
=======================================
  Coverage   88.78%   88.79%           
=======================================
  Files         111      111           
  Lines       17197    17223   +26     
=======================================
+ Hits        15269    15293   +24     
- Misses       1928     1930    +2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Comment on lines -15 to 18
u.Jy/u.sr: u.Unit('W/(m2*sr)'),
u.Jy/u.sr: u.Unit('W/(m2)'),
u.Jy: u.Unit('W/m2'),
u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2*sr)'),
u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2)'),
u.Unit('W/m2/m'): u.Unit('W/m2')
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Hopefully these changes make sense in the context of the sum spectrum using units without per steradian.

Copy link
Contributor

@cshanahan1 cshanahan1 left a comment

Choose a reason for hiding this comment

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

This looks good to me, other than my previous comment about this being specific to per steradians rather than per area, but I think that might be out of scope anyway.

@@ -283,6 +285,73 @@ def standardize_metadata(metadata):
return out_meta


def flux_conversion(spec, values, original_units, target_units):
Copy link
Member

Choose a reason for hiding this comment

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

still don't love this name, but I don't exactly have better suggestions that aren't quite verbose. Since this is internal, I think its fine for now and we can always change it later. Can you maybe add a docstring or longer comment explaining what it does though so future us don't get confused?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Added, let me know if it looks good. If so, I'll merge once CI passes.

Copy link
Member

Choose a reason for hiding this comment

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

thanks!

@javerbukh javerbukh modified the milestones: 3.11, 4.0 Jun 25, 2024
@javerbukh javerbukh merged commit 358d656 into spacetelescope:main Jun 26, 2024
19 checks passed
@javerbukh javerbukh deleted the spec-extract-sum-units branch June 26, 2024 13:00
# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
cubeviz plugin Label for plugins common to multiple configurations testing
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants