Add conical aperture support in the spectral extraction plugin

CHANGES.rst

@@ -13,6 +13,8 @@ New Features
 
 - Live-preview of aperture selection in plugins. [#2664, #2684]
 
+- Add conical aperture support to cubeviz in the spectral extraction plugin [#2679]
+
 Cubeviz
 ^^^^^^^
 

docs/cubeviz/plugins.rst

@@ -291,6 +291,23 @@ Click :guilabel:`EXTRACT` to produce a new 1D spectrum dataset
 from the spectral cube, which has uncertainties propagated by
 `astropy.nddata <https://docs.astropy.org/en/stable/nddata/nddata.html>`_.
 
+If using a simple subset for the spatial aperture, an option to
+make the aperture wavelength dependent will appear. If checked, this will
+create a cone aperture that increases linearly with wavelength.
+The formula for that is::
+
+    radius = ((slice_wavelength / reference_wavelength) *
+        aperture.selected_spatial_region.radius)
+
+The reference wavelength for the cone can be changed using the
+:guilabel:`Adopt Current Slice` button.
+
+The method of aperture masking can also be changed using the
+:guilabel:`Aperture masking method` dropdown. To see a description
+for each of these options, please see
+:ref:`photutils:photutils-aperture-overlap`.
+
+
 .. _cubeviz-aper-phot:
 
 Aperture Photometry

jdaviz/configs/cubeviz/helper.py

@@ -24,6 +24,7 @@ class Cubeviz(ImageConfigHelper, LineListMixin):
 
     _loaded_flux_cube = None
     _loaded_uncert_cube = None
+    _loaded_mask_cube = None
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)

jdaviz/configs/cubeviz/plugins/parsers.py

@@ -254,7 +254,7 @@
 
         if data_type == 'mask':
             # We no longer auto-populate the mask cube into a viewer
-            pass
+            app._jdaviz_helper._loaded_mask_cube = app.data_collection[data_label]
 
         elif data_type == 'uncert':
             app.add_data_to_viewer(uncert_viewer_reference_name, data_label)
@@ -319,6 +319,8 @@
         app._jdaviz_helper._loaded_flux_cube = app.data_collection[data_label]
     elif data_type == 'uncert':
         app._jdaviz_helper._loaded_uncert_cube = app.data_collection[data_label]
+    elif data_type == 'mask':
+        app._jdaviz_helper._loaded_mask_cube = app.data_collection[data_label]
 
 
 def _parse_esa_s3d(app, hdulist, data_label, ext='DATA', flux_viewer_reference_name=None,
@@ -368,8 +370,10 @@
 
     if data_type == 'flux':
         app._jdaviz_helper._loaded_flux_cube = app.data_collection[data_label]
-    if data_type == 'uncert':
+    elif data_type == 'uncert':
         app._jdaviz_helper._loaded_uncert_cube = app.data_collection[data_label]
+    elif data_type == 'mask':
+        app._jdaviz_helper._loaded_mask_cube = app.data_collection[data_label]
 
 
 def _parse_spectrum1d_3d(app, file_obj, data_label=None,
@@ -422,7 +426,9 @@
         elif attr == 'uncertainty':
             app.add_data_to_viewer(uncert_viewer_reference_name, cur_data_label)
             app._jdaviz_helper._loaded_uncert_cube = app.data_collection[cur_data_label]
-        # We no longer auto-populate the mask cube into a viewer
+        elif attr == 'mask':
+            # We no longer auto-populate the mask cube into a viewer
+            app._jdaviz_helper._loaded_mask_cube = app.data_collection[cur_data_label]
 
 
 def _parse_spectrum1d(app, file_obj, data_label=None, spectrum_viewer_reference_name=None):
@@ -470,6 +476,8 @@
     elif data_type == 'uncert':
         app.add_data_to_viewer(uncert_viewer_reference_name, data_label)
         app._jdaviz_helper._loaded_uncert_cube = app.data_collection[data_label]
+    elif data_type == 'mask':
+        app._jdaviz_helper._loaded_mask_cube = app.data_collection[data_label]
 
 
 def _parse_gif(app, file_obj, data_label=None, flux_viewer_reference_name=None,

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

@@ -10,6 +10,8 @@
     NDDataArray, StdDevUncertainty, NDUncertainty
 )
 from traitlets import Any, Bool, Dict, Float, List, Unicode, observe
+from photutils.aperture import CircularAperture
+from specutils import Spectrum1D
 
 from jdaviz.core.custom_traitlets import FloatHandleEmpty
 from jdaviz.core.events import SnackbarMessage, SliceWavelengthUpdatedMessage
@@ -48,14 +50,18 @@
       Subset to use for the spectral extraction, or ``Entire Cube``.
     * ``add_results`` (:class:`~jdaviz.core.template_mixin.AddResults`)
     * :meth:`collapse`
+    * ``wavelength_dependent``:
+      When true, the cone_aperture method will be used to determine the mask.
+    * ``reference_wavelength``:
+      The wavelength that will be used to calculate the radius of the cone through the cube.
+    * ``aperture_method`` (:class:`~jdaviz.core.template_mixin.SelectPluginComponent`):
+      Extract spectrum using an aperture masking method in place of the subset mask.
     """
     template_file = __file__, "spectral_extraction.vue"
     uses_active_status = Bool(True).tag(sync=True)
 
-    # feature flag for cone support
-    dev_cone_support = Bool(False).tag(sync=True)  # when enabling: add entries to docstring
+    # feature flag for background cone support
     dev_bg_support = Bool(False).tag(sync=True)  # when enabling: add entries to docstring
-    dev_subpixel_support = Bool(False).tag(sync=True)  # when enabling: add entries to docstring
 
     active_step = Unicode().tag(sync=True)
 
@@ -69,14 +75,15 @@
     bg_scale_factor = Float(1).tag(sync=True)
     bg_wavelength_dependent = Bool(False).tag(sync=True)
 
-    subpixel = Bool(False).tag(sync=True)
-
     function_items = List().tag(sync=True)
     function_selected = Unicode('Sum').tag(sync=True)
     filename = Unicode().tag(sync=True)
     extracted_spec_available = Bool(False).tag(sync=True)
     overwrite_warn = Bool(False).tag(sync=True)
 
+    aperture_method_items = List(['exact', 'subpixel', 'center']).tag(sync=True)
+    aperture_method_selected = Unicode('exact').tag(sync=True)
+
     # export_enabled controls whether saving to a file is enabled via the UI.  This
     # is a temporary measure to allow server-installations to disable saving server-side until
     # saving client-side is supported
@@ -113,6 +120,12 @@
             selected='function_selected',
             manual_options=['Mean', 'Min', 'Max', 'Sum']
         )
+        self.aperture_method = SelectPluginComponent(
+            self,
+            items='aperture_method_items',
+            selected='aperture_method_selected',
+            manual_options=['exact', 'subpixel', 'center']
+        )
         self._set_default_results_label()
         self.add_results.viewer.filters = ['is_spectrum_viewer']
 
@@ -135,13 +148,11 @@
     @property
     def user_api(self):
         expose = ['function', 'spatial_subset', 'aperture',
-                  'add_results', 'collapse_to_spectrum']
-        if self.dev_cone_support:
-            expose += ['wavelength_dependent', 'reference_wavelength']
+                  'add_results', 'collapse_to_spectrum',
+                  'wavelength_dependent', 'reference_wavelength',
+                  'aperture_method']
         if self.dev_bg_support:
             expose += ['background', 'bg_wavelength_dependent']
-        if self.dev_subpixel_support:
-            expose += ['subpixel']
 
         return PluginUserApi(self, expose=expose)
 
@@ -205,30 +216,43 @@
         """
         spectral_cube = self._app._jdaviz_helper._loaded_flux_cube
         uncert_cube = self._app._jdaviz_helper._loaded_uncert_cube
+        uncertainties = None
 
         # This plugin collapses over the *spatial axes* (optionally over a spatial subset,
         # defaults to ``No Subset``). Since the Cubeviz parser puts the fluxes
         # and uncertainties in different glue Data objects, we translate the spectral
         # cube and its uncertainties into separate NDDataArrays, then combine them:
-        if self.aperture.selected != self.aperture.default_text:
+        if self.aperture.selected != self.aperture.default_text and self.wavelength_dependent:
             nddata = spectral_cube.get_subset_object(
                 subset_id=self.aperture.selected, cls=NDDataArray
             )
-            uncertainties = uncert_cube.get_subset_object(
-                subset_id=self.aperture.selected, cls=StdDevUncertainty
+            if uncert_cube:
+                uncertainties = uncert_cube.get_subset_object(
+                    subset_id=self.aperture.selected, cls=StdDevUncertainty
+                )
+            mask = self.cone_aperture()
+
+        elif self.aperture.selected != self.aperture.default_text:
+            nddata = spectral_cube.get_subset_object(
+                subset_id=self.aperture.selected, cls=NDDataArray
             )
+            if uncert_cube:
+                uncertainties = uncert_cube.get_subset_object(
+                    subset_id=self.aperture.selected, cls=StdDevUncertainty
+                )
+            mask = nddata.mask
         else:
             nddata = spectral_cube.get_object(cls=NDDataArray)
-            uncertainties = uncert_cube.get_object(cls=StdDevUncertainty)
-
+            if uncert_cube:
+                uncertainties = uncert_cube.get_object(cls=StdDevUncertainty)
+            mask = nddata.mask
         # Use the spectral coordinate from the WCS:
         if '_orig_spec' in spectral_cube.meta:
             wcs = spectral_cube.meta['_orig_spec'].wcs.spectral
         else:
             wcs = spectral_cube.coords.spectral
 
         flux = nddata.data << nddata.unit
-        mask = nddata.mask
 
         nddata_reshaped = NDDataArray(
             flux, mask=mask, uncertainty=uncertainties, wcs=wcs, meta=nddata.meta
@@ -284,6 +308,38 @@
 
         return collapsed_spec
 
+    def cone_aperture(self):
+        if not self._app._jdaviz_helper._loaded_mask_cube:
+            snackbar_message = SnackbarMessage(
+                "Cannot create cone aperture without valid mask cube loaded.",
+                color="error",
+                sender=self)
+            self.hub.broadcast(snackbar_message)
+            return
+        # Retrieve mask cube and create array to represent the cone mask
+        mask_cube = self._app._jdaviz_helper._loaded_mask_cube.get_object(cls=Spectrum1D,
+                                                                          statistic=None)
+        flux_cube = self._app._jdaviz_helper._loaded_flux_cube.get_object(cls=Spectrum1D,
+                                                                          statistic=None)
+        masks_boolean_values = np.zeros_like(mask_cube.flux.value)
+
+        # Center is reverse coordinates
+        center = (self.aperture.selected_spatial_region.center.y,
+                  self.aperture.selected_spatial_region.center.x)
+        # TODO: Replace with code for retrieving display_unit in cubeviz when it is available
+        display_unit = flux_cube.spectral_axis.unit
+
+        # Loop through cube and create cone aperture at each wavelength. Then convert that to a
+        # mask using the selected aperture method and add that to a mask cube.
+        for index, wavelength in enumerate(mask_cube.spectral_axis):
+            radius = ((wavelength.to(display_unit).value / self.reference_wavelength) *
+                      self.aperture.selected_spatial_region.radius)
+            aperture = CircularAperture(center, r=radius)
+            slice_mask = aperture.to_mask(method=self.aperture_method_selected).to_image(
+                (len(mask_cube.flux), len(mask_cube.flux[0])))
+            masks_boolean_values[:, :, index] = ~(slice_mask > 0)
+        return masks_boolean_values
+
     def vue_spectral_extraction(self, *args, **kwargs):
         self.collapse_to_spectrum(add_data=True)
 

jdaviz/configs/cubeviz/plugins/spectral_extraction/spectral_extraction.vue

@@ -19,13 +19,13 @@
         hint="Select a spatial region to extract its spectrum."
       />
 
-      <v-row v-if="aperture_selected !== 'Entire Cube' && !aperture_selected_validity.is_aperture && dev_cone_support">
+      <v-row v-if="aperture_selected !== 'Entire Cube' && !aperture_selected_validity.is_aperture">
         <span class="v-messages v-messages__message text--secondary">
             {{aperture_selected}} does not support wavelength dependence (cone support): {{aperture_selected_validity.aperture_message}}.
         </span>
       </v-row>
 
-      <div v-if="aperture_selected_validity.is_aperture && dev_cone_support">
+      <div v-if="aperture_selected_validity.is_aperture">
         <v-row>
           <v-switch
             v-model="wavelength_dependent"
@@ -56,7 +56,7 @@
           <v-row justify="end">
             <j-tooltip tooltipcontent="Select the slice nearest the reference wavelength">
               <plugin-action-button :results_isolated_to_plugin="true" @click="goto_reference_wavelength">
-                Slice to Reference Wavelength
+                Slice to Wavelength
               </plugin-action-button>
             </j-tooltip>
           </v-row>
@@ -80,16 +80,15 @@
         </span>
       </v-row>
 
-      <v-row v-if="bg_selected !== 'None' && !bg_selected_validity.is_aperture && dev_cone_support">
+      <v-row v-if="bg_selected !== 'None' && !bg_selected_validity.is_aperture">
         <span class="v-messages v-messages__message text--secondary">
             {{bg_selected}} does not support wavelength dependence (cone support): {{bg_selected_validity.aperture_message}}.
         </span>
       </v-row>
 
       <div v-if="aperture_selected_validity.is_aperture
                  && bg_selected_validity.is_aperture
-                 && wavelength_dependent
-                 && dev_cone_support">
+                 && wavelength_dependent">
         <v-row>
           <v-switch
             v-model="bg_wavelength_dependent"
@@ -118,28 +117,30 @@
     <div @mouseover="() => active_step='ext'">
       <j-plugin-section-header :active="active_step==='ext'">Extract</j-plugin-section-header>
 
-      <v-row v-if="aperture_selected !== 'None' && !aperture_selected_validity.is_aperture && dev_subpixel_support">
+      <v-row v-if="aperture_selected !== 'None' && !aperture_selected_validity.is_aperture">
         <span class="v-messages v-messages__message text--secondary">
             Aperture: {{aperture_selected}} does not support subpixel: {{aperture_selected_validity.aperture_message}}.
         </span>
       </v-row>
-      <v-row v-if="bg_selected !== 'None' && !bg_selected_validity.is_aperture && dev_subpixel_support">
+      <v-row v-if="bg_selected !== 'None' && !bg_selected_validity.is_aperture">
         <span class="v-messages v-messages__message text--secondary">
             Background: {{bg_selected}} does not support subpixel: {{bg_selected_validity.aperture_message}}.
         </span>
       </v-row>
 
 
       <div v-if="(aperture_selected === 'Entire Cube' || aperture_selected_validity.is_aperture)
-                 && (bg_selected === 'None' || bg_selected_validity.is_aperture)
-                 && dev_subpixel_support">
+                 && (bg_selected === 'None' || bg_selected_validity.is_aperture)">
         <v-row>
-          <v-switch
-            v-model="subpixel"
-            label="Subpixel"
-            hint="Extract spectrum using a subpixel aperture in place of the subset mask"
-            persistent-hint>
-          </v-switch>
+          <v-select
+            :menu-props="{ left: true }"
+            attach
+            :items="aperture_method_items.map(i => i.label)"
+            v-model="aperture_method_selected"
+            label="Aperture masking method"
+            hint="Extract spectrum using an aperture masking method in place of the subset mask."
+            persistent-hint
+            ></v-select>
         </v-row>
       </div>
 

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

@@ -231,3 +231,35 @@ def test_aperture_markers(cubeviz_helper, spectrum1d_cube):
         extract_plg._obj.vue_adopt_slice_as_reference()
         extract_plg._obj.vue_goto_reference_wavelength()
         assert slice_plg.slice == 1
+
+
+def test_cone_aperture(cubeviz_helper, spectrum1d_cube_largest):
+    cubeviz_helper.load_data(spectrum1d_cube_largest)
+    cubeviz_helper.load_regions([CirclePixelRegion(PixCoord(1, 1), radius=0.5)])
+
+    mask_cube = Spectrum1D(flux=np.ones_like(spectrum1d_cube_largest.flux),
+                           spectral_axis=spectrum1d_cube_largest.spectral_axis)
+    cubeviz_helper.load_data(mask_cube, override_cube_limit=True)
+    cubeviz_helper._loaded_mask_cube = cubeviz_helper.app.data_collection[-1]
+
+    extract_plg = cubeviz_helper.plugins['Spectral Extraction']
+    slice_plg = cubeviz_helper.plugins['Slice']
+
+    extract_plg.aperture = 'Subset 1'
+    extract_plg.wavelength_dependent = True
+    assert cubeviz_helper._loaded_mask_cube.get_object(cls=Spectrum1D, statistic=None)
+
+    slice_plg.slice = 1
+    extract_plg._obj.vue_adopt_slice_as_reference()
+    cone_aperture = extract_plg._obj.cone_aperture()
+    assert (cone_aperture.shape ==
+            cubeviz_helper._loaded_flux_cube.get_object(cls=Spectrum1D, statistic=None).shape)
+
+    # Make sure that the cone created when the reference slice is 988 is different
+    # to the cone made at reference slice 1.
+    slice_plg.slice = 988
+    extract_plg._obj.vue_adopt_slice_as_reference()
+    cone_aperture_2 = extract_plg._obj.cone_aperture()
+
+    with pytest.raises(AssertionError, match="Arrays are not equal"):
+        assert np.testing.assert_array_equal(cone_aperture, cone_aperture_2)

jdaviz/conftest.py

@@ -238,6 +238,11 @@ def spectrum1d_cube_larger():
     return _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy, shape=(SPECTRUM_SIZE, 2, 4))
 
 
+@pytest.fixture
+def spectrum1d_cube_largest():
+    return _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy, shape=(1000, 2, 4))
+
+
 @pytest.fixture
 def spectrum1d_cube_custom_fluxunit():
     return _create_spectrum1d_cube_with_fluxunit