Perform corrections to Line Analysis Flux Results

CHANGES.rst

@@ -6,6 +6,9 @@ New Features
 
 - Profile viewers now support plotting with profiles "as steps". [#1595]
 
+- Line flux in the Line Analysis plugin are reported in W/m2 if Spectral Flux is given
+  in Jy [#1564]
+
 Cubeviz
 ^^^^^^^
 

docs/specviz/plugins.rst

@@ -234,6 +234,8 @@ secondary region can be created and selected as the region to fit the linear con
 The statistics returned include the line centroid, gaussian sigma width, gaussian FWHM,
 total flux, and equivalent width.
 
+The line flux results are automatically converted to Watts/meter^2, when appropriate.
+
 Redshift from Centroid
 ----------------------
 

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

@@ -7,7 +7,7 @@
 from glue_jupyter.common.toolbar_vuetify import read_icon
 from traitlets import Bool, List, Float, Unicode, observe
 from astropy import units as u
-from specutils import analysis
+from specutils import analysis, Spectrum1D
 from specutils.manipulation import extract_region
 
 from jdaviz.core.custom_traitlets import FloatHandleEmpty
@@ -332,7 +332,56 @@ def _uncertainty(result):
         for function in FUNCTIONS:
             # TODO: update specutils to allow ALL analysis to take regions and continuum so we
             # don't need these if statements
-            if function == "Equivalent Width":
+            if function == "Line Flux":
+                flux_unit = spec_subtracted.flux.unit
+                # If the flux unit is equivalent to Jy, or Jy per spaxel for Cubeviz,
+                # enforce integration in frequency space
+                if (flux_unit.is_equivalent(u.Jy) or
+                        flux_unit.is_equivalent(u.Jy/u.sr)):
+                    # Perform integration in frequency space
+                    freq_spec = Spectrum1D(
+                        spectral_axis=spec_subtracted.spectral_axis.to(u.Hz,
+                                                                       equivalencies=u.spectral()),
+                        flux=spec_subtracted.flux,
+                        uncertainty=spec_subtracted.uncertainty)
+
+                    raw_result = analysis.line_flux(freq_spec)
+                    # When flux is equivalent to Jy, lineflux result should be shown in W/m2
+                    if flux_unit.is_equivalent(u.Jy/u.sr):
+                        final_unit = u.Unit('W/m2/sr')
+                    else:
+                        final_unit = u.Unit('W/m2')
+
+                    temp_result = raw_result.to(final_unit)
+                    if getattr(raw_result, 'uncertainty', None) is not None:
+                        temp_result.uncertainty = raw_result.uncertainty.to(final_unit)
+
+                # If the flux unit is instead equivalent to power density
+                # (Jy, but defined in wavelength), enforce integration in wavelength space
+                elif (flux_unit.is_equivalent(u.Unit('W/m2/m')) or
+                        flux_unit.is_equivalent(u.Unit('W/m2/m/sr'))):
+                    # Perform integration in wavelength space using MKS unit (meters)
+                    wave_spec = Spectrum1D(
+                        spectral_axis=spec_subtracted.spectral_axis.to(u.m,
+                                                                       equivalencies=u.spectral()),
+                        flux=spec_subtracted.flux,
+                        uncertainty=spec_subtracted.uncertainty)
+                    raw_result = analysis.line_flux(wave_spec)
+                    # When flux is equivalent to Jy, lineflux result should be shown in W/m2
+                    if flux_unit.is_equivalent(u.Unit('W/m2/m'/u.sr)):
+                        final_unit = u.Unit('W/m2/sr')
+                    else:
+                        final_unit = u.Unit('W/m2')
+
+                    temp_result = raw_result.to(final_unit)
+                    if getattr(raw_result, 'uncertainty', None) is not None:
+                        temp_result.uncertainty = raw_result.uncertainty.to(final_unit)
+
+                # Otherwise, just rely on the default specutils line_flux result
+                else:
+                    temp_result = analysis.line_flux(spec_subtracted)
+
+            elif function == "Equivalent Width":
                 if np.any(continuum <= 0):
                     temp_results.append({'function': function,
                                          'result': '',

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

@@ -136,7 +136,7 @@ def test_continuum_surrounding_spectral_subset(specviz_helper, spectrum1d):
     plugin.width = 3
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), 350.89288537581467, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), 2.153181e-13, atol=1e-15)
 
 
 def test_continuum_spectral_same_value(specviz_helper, spectrum1d):
@@ -215,7 +215,7 @@ def test_continuum_subset_spectral_entire(specviz_helper, spectrum1d):
     plugin.width = 3
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), -467.6854635447396, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), -2.79572e-13, atol=1e-15)
 
 
 def test_continuum_subset_spectral_subset2(specviz_helper, spectrum1d):
@@ -248,7 +248,7 @@ def test_continuum_subset_spectral_subset2(specviz_helper, spectrum1d):
     plugin.width = 3
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), 32.520521, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), 1.482418e-14, atol=1e-16)
 
 
 def test_continuum_surrounding_no_right(specviz_helper, spectrum1d):
@@ -276,7 +276,7 @@ def test_continuum_surrounding_no_right(specviz_helper, spectrum1d):
     plugin.width = 3
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), 39.76685499263615, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), 4.204513e-14, atol=1e-16)
 
 
 def test_continuum_surrounding_no_left(specviz_helper, spectrum1d):
@@ -304,7 +304,7 @@ def test_continuum_surrounding_no_left(specviz_helper, spectrum1d):
     plugin.width = 3
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), 146.67186446784513, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), 7.570859e-14, atol=1e-16)
 
 
 def test_subset_changed(specviz_helper, spectrum1d):
@@ -335,4 +335,4 @@ def test_subset_changed(specviz_helper, spectrum1d):
     specviz_helper.app.state.drawer = True
 
     # Values have not yet been validated
-    np.testing.assert_allclose(float(plugin.results[0]['result']), 350.89288537581467, atol=1e-5)
+    np.testing.assert_allclose(float(plugin.results[0]['result']), 2.153181e-13, atol=1e-15)

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

@@ -0,0 +1,123 @@
+from astropy.tests.helper import assert_quantity_allclose
+import astropy.units as u
+from glue.viewers.profile.state import FUNCTIONS as COLLAPSE_FUNCTIONS
+import numpy as np
+import pytest
+from specutils import Spectrum1D
+
+from jdaviz import Cubeviz
+
+import warnings
+warnings.filterwarnings('ignore')
+
+# Maps input spectrum flux unit to expected line analysis flux unit
+expected_lineflux_results = {
+    u.Jy/u.sr: u.Unit('W/(m2*sr)'),
+    u.Jy: u.Unit('W/m2'),
+    u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2*sr)'),
+    u.Unit('W/m2/m'): u.Unit('W/m2')
+}
+
+test_cases = list(expected_lineflux_results.keys())
+
+
+# Test cases for unit gaussian tests
+def _gauss_with_unity_area(x, mean, sigma):
+    ''' Gaussian function with area = 1 unit '''
+    dx = x - mean
+    n = 1/(sigma*np.sqrt(2*np.pi))
+    g = n * np.exp(-dx**2/(2*sigma**2))
+    return g
+
+
+mn = 1.5
+sig = 0.01
+
+# 5 cases: replace fl_wave with fnu_freq, fnu_wave, flam_wave, flam_freq, or fl_wave
+# to try the different cases here
+unit_flux_gaussian_test_cases = []
+
+# 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.0E26*u.Jy
+fnu_freq = Spectrum1D(spectral_axis=freq, flux=flux_freq)
+unit_flux_gaussian_test_cases.append(fnu_freq)
+fnu_wave = Spectrum1D(spectral_axis=fnu_freq.wavelength, flux=flux_freq)
+unit_flux_gaussian_test_cases.append(fnu_wave)
+
+# 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.0*u.W/u.m**2/u.m
+flam_wave = Spectrum1D(spectral_axis=lam, flux=flux_wave)
+unit_flux_gaussian_test_cases.append(flam_wave)
+flam_freq = Spectrum1D(spectral_axis=flam_wave.frequency, flux=flux_wave)
+unit_flux_gaussian_test_cases.append(flam_freq)
+
+
+def _calculate_line_flux(viz_helper):
+    '''
+    Returns the line flux calculation by opening the Line Analysis Plugin
+    Assumes the plugin hasn't been opened yet
+    '''
+    # Open the plugin and force the calculation
+    viz_helper.app.state.drawer = True
+    line_analysis_plugin = viz_helper.app.get_tray_item_from_name('specviz-line-analysis')
+    line_analysis_plugin.open_in_tray()
+    # Retrieve Results
+    for result in line_analysis_plugin.results:
+        if result['function'] == 'Line Flux':
+            return result
+
+
+@pytest.mark.filterwarnings('ignore')
+@pytest.mark.parametrize('spectra_fluxunit', test_cases)
+def test_cubeviz_collapse_fluxunits(spectrum1d_cube_custom_fluxunit, spectra_fluxunit):
+    ''' Calculates line flux and checks the units for each collapse function '''
+    data = spectrum1d_cube_custom_fluxunit(spectra_fluxunit)
+    for function in COLLAPSE_FUNCTIONS:
+        # Initialize Cubeviz with specific data and collapse function
+        cubeviz_helper = Cubeviz()
+        data_label = "Test Cube"
+        cubeviz_helper.load_data(data, data_label=data_label)
+        cubeviz_helper.app.get_viewer('spectrum-viewer').state.function = function
+
+        lineflux_result = _calculate_line_flux(cubeviz_helper)
+
+        autocollapsed_spectrum_unit = cubeviz_helper.specviz.get_spectra()[data_label + "[FLUX]"
+                                                                           ].flux.unit
+        # Futureproofing: Eventually Cubeviz autocollapse will change the flux units of the
+        # spectra depending on whether the spectrum was collapsed OVER SPAXELS or not. Only
+        # do the assertion check if we KNOW what the expected lineflux results should be
+        if autocollapsed_spectrum_unit in expected_lineflux_results.keys():
+            assert u.Unit(lineflux_result['unit']) == expected_lineflux_results[spectra_fluxunit]
+
+
+@pytest.mark.filterwarnings('ignore')
+@pytest.mark.parametrize('test_case', unit_flux_gaussian_test_cases)
+def test_unit_gaussian(specviz_helper, test_case):
+    '''
+    Test an Area 1 Gaussian and ensure the result returns in W/m2
+    Test provided by Patrick Ogle
+    '''
+    specviz_helper.load_data(test_case)
+
+    lineflux_result = _calculate_line_flux(specviz_helper)
+    assert_quantity_allclose(float(lineflux_result['result']) * u.Unit(lineflux_result['unit']),
+                             1*u.Unit('W/m2'))
+
+
+@pytest.mark.filterwarnings('ignore')
+def test_unit_gaussian_mixed_units_per_steradian(specviz_helper):
+    '''
+    A special unit test of Area 1 with mixed units. Should return W/m2sr
+    Test provided by Patrick Ogle
+    '''
+    # 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)*1E3*u.erg/u.s/u.cm**2/u.Angstrom/u.sr
+    fl_wave = Spectrum1D(spectral_axis=lam_a, flux=flx_wave)
+
+    specviz_helper.load_data(fl_wave)
+    lineflux_result = _calculate_line_flux(specviz_helper)
+    assert_quantity_allclose(float(lineflux_result['result']) * u.Unit(lineflux_result['unit']),
+                             1*u.Unit('W/(m2sr)'))

jdaviz/conftest.py

@@ -156,9 +156,9 @@ def spectrum_collection(spectrum1d):
     return result
 
 
-@pytest.fixture
-def spectrum1d_cube():
-    flux = np.arange(16).reshape((2, 2, 4)) * u.Jy
+def _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy):
+
+    flux = np.arange(16).reshape((2, 2, 4)) * fluxunit
     wcs_dict = {"CTYPE1": "RA---TAN", "CTYPE2": "DEC--TAN", "CTYPE3": "WAVE-LOG",
                 "CRVAL1": 205, "CRVAL2": 27, "CRVAL3": 4.622e-7,
                 "CDELT1": -0.0001, "CDELT2": 0.0001, "CDELT3": 8e-11,
@@ -168,6 +168,16 @@ def spectrum1d_cube():
     return Spectrum1D(flux=flux, wcs=w)
 
 
+@pytest.fixture
+def spectrum1d_cube():
+    return _create_spectrum1d_cube_with_fluxunit(u.Jy)
+
+
+@pytest.fixture
+def spectrum1d_cube_custom_fluxunit():
+    return _create_spectrum1d_cube_with_fluxunit
+
+
 @pytest.fixture
 def mos_spectrum1d():
     '''