Bugfix for flux unit conversion in spectral cubes

[bmorris3]

NOTE

This PR has a significant update below. I'm keeping the original description here for traceability.

Original description

Fixes #3086

In (somewhat) recent changes, the default y-axis on the spectrum-viewer in Cubeviz is MJy and the spectral extraction plugin is used directly to produce the default spectrum in the spectrum viewer. If a user makes a spatial subset in Cubeviz that appears to cover every data spaxel, one might expect the live-extracted spectrum to match the default spectral extraction, but it doesn't, see #3086.

This appears to come from the calculation of the pixel area in steradians to convert the flux unit from MJy/sr to MJy. The spectral extraction plugin currently counts every element in the spectral cube towards the sky area, so that the conversion looks like:

(spectral flux density [MJy/sr]) * (full cube pixel area coverage [sr]))

This would be true if every element of the spectral cube contained data. But in general, they don't.

Spectral cubes from JWST/MIRI contain many spaxels near the borders that contain no data at any wavelength, and the sky footprint of the MIRI cube changes from wavelength to wavelength within one spectral cube. Here's Cubeviz with the DQ array visible to show the pixels that do not contain observations:

Those red-colored pixels should not count towards the pixel area, but on main right now, they are included in the conversion from MJy/sr -> MJy. The bug is exposed in #3086 because @rosteen drew a subset that included some but not all non-science pixels. This is a simple demonstration that the approach counting non-science pixels in the pixel area gives you the wrong answer.

This PR excludes "non-science" pixels from the pixel area calculation during flux unit conversion. This brings the default extraction into agreement with the all-pixel-subset extraction:

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.

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)? 🐱

[kecnry]

would there be any use for this logic within aperture.get_mask instead of a property in the plugin?

[bmorris3]

Since ApertureSubsetSelect is meant to be general, I'd advocate for this particular solution to be associated with spectral extraction. There are different ways that you might want to handle NaN masks, for example, in spectral cubes vs images. If you had a flat-fielded image that produced a NaN, you might want your aperture photometry to do something like zero-filling or interpolating over the NaN, in which case you would still count the sky area in that pixel.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 88.78%. Comparing base (73e54d0) to head (c1ae0f2).
Report is 162 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3088      +/-   ##
==========================================
+ Coverage   88.76%   88.78%   +0.01%     
==========================================
  Files         111      111              
  Lines       17246    17276      +30     
==========================================
+ Hits        15308    15338      +30     
  Misses       1938     1938              

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

[rosteen]

This is definitely already an improvement - if I select a rectangular subset that covers the whole cube, I get basically the same spectrum out as the "whole cube" extraction (there are very minor differences) as expected. However, I'm still skeptical of the extraction for a circular subset. In the video below you can see that the circular subset covers basically all of the regions with flux from the source, but the extracted spectrum is about half of the "whole cube" spectrum. I would expect this extraction to have 90%+ of the total flux in it.

Screen.Recording.2024-07-15.at.5.00.52.PM.mov

[bmorris3]

@rosteen - I'm not sure the result is as unexpected as you think. The y-axis in your screengrab doesn't go down to zero, so you aren't seeing half the flux in that extraction. You may be seeing less flux than you would expect if the source were much much brighter than the background, but the background is nonzero. Since the amount of area in the non-selected spaxels is still pretty significant, I'm not sure this is wrong.

[rosteen]

@rosteen - I'm not sure the result is as unexpected as you think. The y-axis in your screengrab doesn't go down to zero, so you aren't seeing half the flux in that extraction.

True, but hovering over the values/looking at the axes, regardless of the zoom it seems like the extracted subset spectrum is ~65% of the full cube spectrum, which is much lower than I'd expect. Looking at a similar subset on 3.10.2, it looks like the subset spectrum is ~85% of the full cube:

I'm open to being wrong about this, but I'm still not convinced. I could take a more rigorous look defining exactly the same subset and dividing the actual Spectrum1D objects if you would like (although we're about to have our hack day and I'm off tomorrow, so it might be delayed).

[rosteen]

Update: using the Subset Tools plugin to define the exact same circular subset (xcen=23.329, ycen=21.471, radius=15.375) on the example notebook data, in 3.10.2 the extracted sum from the subset is (on average) 84.7% of the flux of the full cube extraction, while with this PR it's 55.6%. That seems like a big difference to me.

[bmorris3]

@rosteen I'm also open to being confused, but here's why I think it's right:

• on this branch, the sum of the flux is multiplied by the area of the subset, rather than the area of the full cube. We expect to extract less flux, but how much?
• the ratio of the extracted flux on this branch to the flux on v3.10.2 should be the ratio of the area of science spaxels within the subset (new normalization) to the area of all science spaxels (old normalization)
• that ratio is:

>>> import numpy as np

>>> subset_area = np.pi * 15.375 ** 2
>>> non_nan_area = np.count_nonzero(~np.isnan(data.get_component('flux').data[..., 0]))

>>> subset_area / non_nan_area
0.6702554610807759

You observed in the previous comment that the ratio of flux between this branch and v3.10.2 is 55.6/84.7 = 65.6%. Does that make sense?

[pllim]

Is this affected by the bug I fixed in #2936 ?

[bmorris3]

@pllim – I don't think they're related, but let me know if I'm missing something.

[bmorris3]

After a very helpful discussion led by @rosteen, we've found further bugs to fix in our spectral extraction approach. I've converted this PR back to draft while I work on the rest of the fixes.

[bmorris3]

Significant update

Rather than updating the description above, I'm putting this update here so the history can be rediscovered in the future

@rosteen pointed out in tag up that the approach on main, and tweaked in this PR, does the following conversion to convert the extracted spectrum within an aperture to flux:

sum( (flux [MJy / sr]) * (# pixels) * (solid angle in one pixel [sr / pix]) )

Since the sum is already adding the flux from each pixel, the multiplication by the number of pixels is not necessary. As a result, all extracted spectra have greater fluxes than they should. This was not covered by the tests, which did not yet contain a scientific validation test to compare the extracted flux with externally validated measurements.

The most recent updates to this PR:

• remove the double-counted number of pixels from the integral of the flux within the aperture
• fix the normalization for background extractions
• add scientific validation tests – more on this below

Scientific validation

The new test compares the extracted spectrum from MIRI CH1 IFU cubes for delta UMi (A1Van star) and HD 159222 (G1V star) with CALSPEC model spectra. These model spectra are already flux calibrated using other observatories, and are used to compute the flux calibration for JWST. Current CALSPEC models are available here.

We can compare the spectral extraction results from Cubeviz with these model spectra, without further normalization of the models or data. When comparing any flux-calibrated observation against any model, the expected mismatch can be as large as ~10%. The agreement is much better for stars used to compute the flux calibration, which are the targets selected for this test.

One target gets most of the weight in the flux-calibration used by the jwst pipeline, which is the A1Van star delta UMi. Here's a comparison between the flux-calibrated model spectrum and the extracted spectrum in Cubeviz with the updates in this PR:

Both observations in this test are from MIRI. The test ensures MAD relative flux tolerance better than 10 ppm for delta UMi, and 1% for HD 159222. This could be improved by including fringing corrections in the future – read more on jdox. It's likely that the accuracy of the flux calibration of the data products available on MAST will evolve with time. For the latest updates on MIRI flux calibration, see jdox.

Big thanks to @drlaw1558 and @skendrew for teaching me all about MIRI flux calibration!

[rosteen]

This looks good to me now - thanks for the extensive comments, for doing the science validation, and for the extensive test coverage!

[kecnry]

is aperture_area_along_spectral and bg_area_along_spectral used anywhere or should they be removed?

[bmorris3]

It's needed for scaling the background flux by the ratio of the background to aperture areas:

jdaviz/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py

Line 568 in c1ae0f2

bg_spec *= self.aperture_area_along_spectral / self.bg_area_along_spectral

The last important place to check is this though:

jdaviz/jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.py

Line 516 in c1ae0f2

pix_scale_factor = self.aperture_area_along_spectral * self.spectral_cube.meta.get('PIXAR_SR', 1.0) # noqa

@gibsongreen – do you think we need that aperture scale factor?

[bmorris3]

Confirmed with @gibsongreen that we should leave that term there.

[kecnry]

Just one small question if we can remove unneeded code, otherwise the changes look reasonable to me. Thanks to @rosteen for noticing something was 🐟-y!

[bmorris3]

@rosteen == 🦸🏻‍♂️.