Perform corrections to Line Analysis Flux Results

[duytnguyendtn]

Description

This PR performs corrections to the Line Analysis as requested by @PatrickOgle to get the base units of the Line Flux to Watts/meter^2.

One remaining task is to multiply with the pixel area if provided. I'm not sure how that is represented in the header, so I feel like more guidance is needed to get that working. Even better if an example dataset can be provided with PIXAR_SR is set. Some obvious ones to start with:

1. Is PIXAR_SR always in units of steradians^2?
2. Is PIXAR_SR only provided in cases the flux is in units of X/sr^2? (Effectively allowing the two to cancel?)
3. If (1) is true and (2) is false, does that mean a unit of steradians^2 should be ADDED to the line flux unit reported?

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 change log needed? If yes, is it added to CHANGES.rst?
• Is a milestone set?
• After merge, any internal documentations need updating (e.g., JIRA, Innerspace)?

[codecov]

Codecov Report

Base: 86.28% // Head: 86.39% // Increases project coverage by +0.11% 🎉

Coverage data is based on head (1557d16) compared to base (8050caf).
Patch coverage: 91.30% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1564      +/-   ##
==========================================
+ Coverage   86.28%   86.39%   +0.11%     
==========================================
  Files          94       94              
  Lines        9366     9431      +65     
==========================================
+ Hits         8081     8148      +67     
+ Misses       1285     1283       -2     

Impacted Files	Coverage Δ
...igs/specviz/plugins/line_analysis/line_analysis.py	97.89% <91.30%> (+0.66%)	⬆️
jdaviz/configs/mosviz/plugins/parsers.py	89.54% <0.00%> (+0.08%)	⬆️
jdaviz/configs/imviz/helper.py	96.61% <0.00%> (+0.11%)	⬆️
jdaviz/core/events.py	93.93% <0.00%> (+0.49%)	⬆️
...nfigs/imviz/plugins/links_control/links_control.py	100.00% <0.00%> (+7.14%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[PatrickOgle]

JWST spectra have two relevant keywords:
PIXAR_A2 (arcsecond^2)
PIXAR_SR (steradian) [Not steradian^2, which is not the right units]
The conversion between arcsecond^2 and steradian is sr = 4.255E10 arcsecond^2 (approximately)

JWST spectra will be in either:

1. MJy/sr--- multiplying by PIXAR_SR converts to MJy. This is the correct thing to do
   if the spectrum was summed over spaxels, but not if it is a median or average.
   It is probably not a good idea to guess or try to track this. Rather, there should be a separate PR
   that does the right unit conversion on collapse. This should happen before line analysis.

2. MJy-- Should not be multiplied by solid angle.

So, the bottom line is, do not multiply by pixel solid angle in this PR. This should be handled elsewhere.

[PatrickOgle]

The thing this PR should do is calculate the line flux in units like W/m^2/sr (for inputs like MJy/sr)
or W/m^2 (for inputs like MJy).

[duytnguyendtn]

One small edge case concerning when the Line Analysis units return in MJy * um; I asked for clarification from @PatrickOgle, but otherwise, this should be ready for review! I'll draft up some tests in the meantime

[pllim]

Does user-facing documentation also need updating?

Is this tested anywhere?

[duytnguyendtn]

Is this tested anywhere?

I'll draft up some tests in the meantime

@pllim Yup, tests are incoming!

[pllim]

This PR has a conflict. I think that is why CI is stuck. Please rebase. Thanks!

[kecnry]

Just a few initial comments - I'll take another look once tests are implemented.

If this is the expected behavior, should this logic live in specutils instead of custom jwst logic in jdaviz? Then we can just keep things simple in jdaviz, call specutils directly, and users won't be confused by a discrepancy between jdaviz and native specutils calls? 🤷

[kecnry]

but this else also catches cases where the viewer collapse function is something other than min/max/sum, right? Is that expected behavior for the user or would that be confusing?

[duytnguyendtn]

The path here is super confusing (see #1564 (comment)), so let me try to explain it and let's see if it makes sense, or if I confused myself 😅

If the Flux units are in surface brightness (e.g. has the steradians in it), but PIXAR_SR was not provided, the units should ALWAYS be in W/(m2 sr). It doesn't matter what the collapse function is. The conversion to W/m2 ONLY happens when PIXAR_SR exists, AND the function is either min, max, or sum.

To directly answer your question, if the function is mean/median and the flux is in surface brightness, then the units SHOULD NOT be converted to W/m2, instead SHOULD keep the steradian. Likewise, if the collapse function is min/max/sum and PIXAR_SR was NOT provided, then likewise the units should keep the steradian.

How'd I do?

[duytnguyendtn]

Converting to draft while I handle test failures (my local testing env broke :/)

[duytnguyendtn]

Ok, new update/scope change from @PatrickOgle. I'll try to summarize here, but please correct me if I get something wrong!

So the root of the min/max/sum multiplying by PIXAR_SR is due to the fact that Cubeviz' autocollapse doesn't respect the units properly. Currently, Cubeviz collapses all retain the same unit of the flux, regardless of which function you choose. However, this is technically incorrect. When you integrate OVER SPAXELS (aka min/max/sum), you should introduce a spaxel unit of steradians into the flux, which EFFECTIVELY drops the 1/sr from the MJy/sr. When line analysis computes the flux, it results in the expected behavior above.

What this means is the logic in

jdaviz/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py

Lines 354 to 363 in c2a1126

	# Multiply by PIXAR_SR if available
	viewer = self.app.get_viewer('spectrum-viewer')
	pixar_sr = meta.get('PIXAR_SR',
	meta.get('_primary_header', {}).get('PIXAR_SR', ''))
	if (hasattr(viewer.state, 'function') and pixar_sr and
	# We can rely on the viewer's function because Cubeviz only allows
	# one cube per instance. If this changes, this loop will need to be
	# modified to look at the specific data entry's collapse function
	viewer.state.function.lower() in ('minimum', 'maximum', 'sum')):
	temp_result = temp_result * (float(pixar_sr) * u.Unit('sr'))

should REALLY be happening TO THE FLUX AT AUTOCOLLAPSE. The plugin should not be doing this fix. @PatrickOgle suggested the plugin should respect the units in the spectrum, even though this is technically a bug, and this fix should move up to the autocollapse in the glue logic. @astrofrog should be made aware of this.

As a result, I've moved this PR back into draft until I can make those changes. Will mark it ready once I make those changes/update tests.

In addition, the integration of Hz is a candidate to be moved up to specutils, but for the time being will remain here. It is possible that the two integrations (over wave vs. over freq) are producing the same results, but we aren't sure of that. A further investigation should be made to figure out if there are any assumptions being made (like if the frequency is in the center of the subset), and verify the results are correct.

[pllim]

I think #1564 (comment) makes sense. But let me make sure I understand: Are you saying that the real fix is in specutils but for now we need to workaround it in our own collapse method? Is there an upstream issue about this in specutils?

[duytnguyendtn]

@pllim There are two "upstream" fixes that should be considered for this PR's next steps:

1. The first concerns the autocollapsed spectrum in Cubeviz, which needs to respect the spaxel units for the integration. This logic was mentioned to live in Glue, hence why I pinged Tom for this one
2. Calculating the Line Flux in frequency space might need to live in specutils depending on the implementation. If calculating the line flux in both frequency and wavelength space result in the same answer, then we're okay. But if there's a different, my understanding is the frequency space approach is preferred. I'm a little shaky on this one, maybe @PatrickOgle can give recount the reasoning here?

[pllim]

Just talking about glue-core , the collapse function is at https://github.com/glue-viz/glue/blob/f8c1e9646e7252aa13f5f58c4838cb9b670875d3/glue/utils/array.py#L431 and looks like data has shed unit at that point. I don't see why glue should care about astronomy-specific computations so deep down the engine. Unless Tom R has a better idea, I think we need to keep track of these units somewhere in the Cubeviz parser logic but I got lost in the code.

This is the only place we are calling compute_statistic directly and seems like it only is for Subset, so maybe not even used in the initial loading:

jdaviz/jdaviz/app.py

Lines 835 to 840 in c2a9b82

	# Compute reduced data based on the current viewer's statistic
	# function. This doesn't seem particularly useful, but better
	# to be consistent.
	reduced_data = Data(x=value.data.compute_statistic(
	stat_func, value.data.id[xid],
	subset_state=value.subset_state.copy(), axis=1))

This triggers the spectrum collapse somehow but where is the actual code? I cannot find it. Does it magically triggers jdaviz/configs/specviz/plugins/parsers.py?

jdaviz/jdaviz/configs/cubeviz/plugins/parsers.py

Lines 192 to 193 in 9d7bc4a

	# Add flux to spectrum viewer
	app.add_data_to_viewer('spectrum-viewer', data_label)

I see this but I think it only targets wavelength unit, not flux:

jdaviz/jdaviz/configs/cubeviz/plugins/parsers.py

Line 91 in 9d7bc4a

def _return_spectrum_with_correct_units(flux, wcs, metadata, data_type, target_wave_unit=None,

[pllim]

I might have tracked down the collapse code to here:

jdaviz/jdaviz/configs/default/plugins/plot_options/plot_options.py

Lines 140 to 141 in 9d7bc4a

	self.collapse_function = PlotOptionsSyncState(self, self.viewer, self.layer, 'function',
	'collapse_func_value', 'collapse_func_sync')

and here:

https://github.com/glue-viz/glue/blob/f8c1e9646e7252aa13f5f58c4838cb9b670875d3/glue/viewers/profile/state.py#L332

Yeah, not sure how to ask it to handle unit there. That code is fully generalized to handle any given attribute in Data.

Maybe all we need to do is to change this line to show correct collapsed unit on the axis:

jdaviz/jdaviz/configs/specviz/plugins/viewers.py

Line 492 in 5d09c46

self.figure.axes[1].label = f"{flux_unit_type} [{data.flux.unit.to_string()}]"

[PatrickOgle]

Here is the full test suite, including the mixed units and 1/sr:

def gauss_with_unity_area(x, mean, sigma):
dx = x - mean
n = 1/(sigmanp.sqrt(2np.pi))
g = n np.exp(-dx**2/(2sigma**2))
return g
mn = 1.5
sig = 0.01

Unit-flux gaussian in frequency space

freq = np.arange(1, 2, 0.001)*u.Hz
flux_freq = gauss_with_unity_area(freq.value, mn, sig)1.0E26u.Jy
fnu_freq = Spectrum1D(spectral_axis=freq, flux=flux_freq)
fnu_wave = Spectrum1D(spectral_axis = fnu_freq.wavelength, flux=flux_freq)

Unit-flux gaussian in wavelength space

lam = np.arange(1,2, 0.001)*u.m
flux_wave = gauss_with_unity_area(lam.value, mn, sig)1.0u.W/u.m**2/u.m
flam_wave = Spectrum1D(spectral_axis=lam, flux=flux_wave)
flam_freq = Spectrum1D(spectral_axis = flam_wave.frequency, flux=flux_wave)

Unit-flux gaussian in wavelength space, mixed units, per steradian

lam_a = np.arange(1,2, 0.001)*u.Angstrom
flx_wave = gauss_with_unity_area(lam_a.value, mn, sig)1E3u.erg/u.s/u.cm**2/u.Angstrom/u.sr
fl_wave = Spectrum1D(spectral_axis=lam_a, flux=flx_wave)

Display Spectrum 1 D in Specviz
(5 cases: replace fl_wave with fnu_freq, fnu_wave, flam_wave, flam_freq, or fl_wave to try the different cases below)

from jdaviz import Specviz
specviz = Specviz()
specviz.load_data(fl_wave)
specviz.app

Then open Line Analysis Plugin to get Flux = 1 W/m^2/[sr] in all 5 cases.

[PatrickOgle]

Happy to give a final review once the above unit tests are in place.

[duytnguyendtn]

@PatrickOgle Ok done!

[PatrickOgle]

Gives the correct answer for unit area Gaussian in wavelength, frequency space, and mixed units.

[kecnry]

Works for me now, thanks (assuming CI will pass)!

Does the "upstream fix required" label still apply? Do we need a ticket to consider whether to move this logic to specutils as a follow-up or did we decide it doesn't belong there at all?

[kecnry]

hopefully this will be improved upstream some day... but does the trick for now!

[duytnguyendtn]

I created astropy/specutils#980 in specutils to make sure we don't forget to discuss

[kecnry]

my comment on these lines was referring to the fact that .uncertainty (or any arbitrary attribute added to a quantity) isn't supported by astropy and so is stripped from the quantity object when calling .to (see astropy/specutils#961 for some related discussion).

[duytnguyendtn]

Ahh right. Sorry my comment was meant as a direct response to your general review. Sorry, poor placement

[duytnguyendtn]

Thanks for the eyes @pllim @kecnry @PatrickOgle! I appreciate everyone's patience as we hammered this PR through its multiple iterations!