Handle aperture photometry when images linked by WCS for advanced cas…

…es (#2154)

* EXP: Build test case for
advanced aperture photometry.

Add tests that will fail without fix

Remove warning from doc

* FEAT: Add parent info to ROI subset

BUG: Make sure recentering works in Subset Tools

* BUG: Fix aperture incorrectly applied
to images linked by WCS if sky projection is different

* Fix PEP 8 warning
and remove outdated comment

* Add tests for sky bg
and remove debug comment

* Apply suggestions from code reviews

---------

Co-authored-by: Brett M. Morris <[email protected]>

CHANGES.rst

@@ -6,6 +6,8 @@ New Features
 
 - Improved design of Launcher and pass filepath arg from cli when no config specified. [#2311]
 
+- Subset Tools plugin now displays the parent data of a spatial (ROI) subset. [#2154]
+
 Cubeviz
 ^^^^^^^
 
@@ -34,6 +36,16 @@ Cubeviz
 Imviz
 ^^^^^
 
+- Fixed Subset Tools unable to re-center non-composite spatial subset on an image
+  that is not the reference data when linked by WCS. [#2154]
+
+- Fixed inaccurate results when aperture photometry is performed on non-reference data
+  that are of a different pixel scale or are rotated w.r.t. the reference data when
+  linked by WCS. [#2154]
+
+- Fixed wrong angle translations between sky regions in ``regions`` and ``photutils``.
+  They were previously off by 90 degrees. [#2154]
+
 Mosviz
 ^^^^^^
 

docs/imviz/plugins.rst

@@ -165,10 +165,6 @@ Simple Aperture Photometry
 
     Regardless of your workflow, any WCS distortion in an image is ignored.
 
-    When you have multiple images linked by WCS and they have different
-    pixel scales or sky orientation, the selected aperture may be incorrectly defined
-    for images that are not the reference image (usually the first loaded one).
-
 This plugin performs simple aperture photometry
 and plots a radial profile for one object within
 an interactively selected region. A typical workflow is as follows:

jdaviz/configs/default/plugins/subset_plugin/subset_plugin.py

@@ -126,23 +126,10 @@ def _sync_selected_from_ui(self, change):
         if m != self.session.edit_subset_mode.edit_subset:
             self.session.edit_subset_mode.edit_subset = m
 
-    '''
-    # This will be needed once we use a dropdown instead of the actual
-    # g-subset-mode component
-    @observe("mode_selected")
-    def _mode_selected_changed(self, event={}):
-        if self.session.edit_subset_mode != self.mode_selected:
-            self.session.edit_subset_mode = self.mode_selected
-    '''
-
     def _unpack_get_subsets_for_ui(self):
         """
         Convert what app.get_subsets returns into something the UI of this plugin
         can display.
-
-        Returns
-        -------
-
         """
         subset_information = self.app.get_subsets(self.subset_selected, simplify_spectral=False,
                                                   use_display_units=True)
@@ -164,12 +151,18 @@ def _unpack_get_subsets_for_ui(self):
             subset_state = spec["subset_state"]
             glue_state = spec["glue_state"]
             if isinstance(subset_state, RoiSubsetState):
+                subset_definition.append({
+                    "name": "Parent", "att": "parent",
+                    "value": subset_state.xatt.parent.label,
+                    "orig": subset_state.xatt.parent.label})
+
                 if isinstance(subset_state.roi, CircularROI):
                     x, y = subset_state.roi.center()
                     r = subset_state.roi.radius
-                    subset_definition = [{"name": "X Center", "att": "xc", "value": x, "orig": x},
-                                         {"name": "Y Center", "att": "yc", "value": y, "orig": y},
-                                         {"name": "Radius", "att": "radius", "value": r, "orig": r}]
+                    subset_definition += [
+                        {"name": "X Center", "att": "xc", "value": x, "orig": x},
+                        {"name": "Y Center", "att": "yc", "value": y, "orig": y},
+                        {"name": "Radius", "att": "radius", "value": r, "orig": r}]
 
                 elif isinstance(subset_state.roi, RectangularROI):
                     for att in ("Xmin", "Xmax", "Ymin", "Ymax"):
@@ -186,7 +179,7 @@ def _unpack_get_subsets_for_ui(self):
                     rx = subset_state.roi.radius_x
                     ry = subset_state.roi.radius_y
                     theta = np.around(np.degrees(subset_state.roi.theta), decimals=_around_decimals)
-                    subset_definition = [
+                    subset_definition += [
                         {"name": "X Center", "att": "xc", "value": xc, "orig": xc},
                         {"name": "Y Center", "att": "yc", "value": yc, "orig": yc},
                         {"name": "X Radius", "att": "radius_x", "value": rx, "orig": rx},
@@ -197,12 +190,12 @@ def _unpack_get_subsets_for_ui(self):
                     x, y = subset_state.roi.center()
                     inner_r = subset_state.roi.inner_radius
                     outer_r = subset_state.roi.outer_radius
-                    subset_definition = [{"name": "X Center", "att": "xc", "value": x, "orig": x},
-                                         {"name": "Y Center", "att": "yc", "value": y, "orig": y},
-                                         {"name": "Inner radius", "att": "inner_radius",
-                                          "value": inner_r, "orig": inner_r},
-                                         {"name": "Outer radius", "att": "outer_radius",
-                                          "value": outer_r, "orig": outer_r}]
+                    subset_definition += [{"name": "X Center", "att": "xc", "value": x, "orig": x},
+                                          {"name": "Y Center", "att": "yc", "value": y, "orig": y},
+                                          {"name": "Inner radius", "att": "inner_radius",
+                                           "value": inner_r, "orig": inner_r},
+                                          {"name": "Outer radius", "att": "outer_radius",
+                                           "value": outer_r, "orig": outer_r}]
 
                 subset_type = subset_state.roi.__class__.__name__
 
@@ -282,6 +275,9 @@ def vue_update_subset(self, *args):
                 return
             sub_states = self.subset_states[index]
             for d_att in sub:
+                if d_att["att"] == 'parent':  # Read-only
+                    continue
+
                 if d_att["att"] == 'theta':  # Humans use degrees but glue uses radians
                     d_val = np.radians(d_att["value"])
                 else:
@@ -361,23 +357,33 @@ def vue_recenter_subset(self, *args):
             raise NotImplementedError(
                 f'Cannot recenter: is_centerable={self.is_centerable}, config={self.config}')
 
+        from astropy.wcs.utils import pixel_to_pixel
         from photutils.aperture import ApertureStats
         from jdaviz.core.region_translators import regions2aperture, _get_region_from_spatial_subset
 
         try:
-            reg = _get_region_from_spatial_subset(self, self.subset_selected)
+            reg = _get_region_from_spatial_subset(self, self.subset_select.selected_subset_state)
             aperture = regions2aperture(reg)
             data = self.dataset.selected_dc_item
             comp = data.get_component(data.main_components[0])
             comp_data = comp.data
-            phot_aperstats = ApertureStats(comp_data, aperture)
-
-            # Centroid was calculated in selected data, which might or might not be
-            # the reference data. However, Subset is always defined w.r.t.
-            # the reference data, so we need to convert back.
-            viewer = self.app._jdaviz_helper.default_viewer
-            x, y, _, _ = viewer._get_real_xy(
-                data, phot_aperstats.xcentroid, phot_aperstats.ycentroid, reverse=True)
+            phot_aperstats = ApertureStats(comp_data, aperture, wcs=data.coords)
+
+            # Sky region from WCS linking, need to convert centroid back to pixels.
+            if hasattr(reg, "to_pixel"):
+                # Centroid was calculated in selected data.
+                # However, Subset is always defined w.r.t. its parent,
+                # so we need to convert back.
+                x, y = pixel_to_pixel(
+                    data.coords,
+                    self.subset_select.selected_subset_state.xatt.parent.coords,
+                    phot_aperstats.xcentroid,
+                    phot_aperstats.ycentroid)
+
+            else:
+                x = phot_aperstats.xcentroid
+                y = phot_aperstats.ycentroid
+
             if not np.all(np.isfinite((x, y))):
                 raise ValueError(f'Invalid centroid ({x}, {y})')
         except Exception as err:

jdaviz/configs/default/plugins/subset_plugin/subset_plugin.vue

@@ -78,7 +78,14 @@
       </v-row>
 
       <v-row v-for="(item, index2) in region" class="row-no-outside-padding">
-        <v-text-field
+        <v-text-field v-if="item.name === 'Parent'"
+          :label="item.name"
+          :value="item.value"
+          style="padding-top: 0px; margin-top: 0px"
+          :readonly="true"
+          hint="Subset was defined with respect to this reference data (read-only)"
+        ></v-text-field>
+        <v-text-field v-if="item.name !== 'Parent'"
           :label="item.name"
           v-model.number="item.value"
           type="number"

jdaviz/configs/imviz/helper.py

@@ -215,6 +215,48 @@ def link_data(self, **kwargs):
         """
         link_image_data(self.app, **kwargs)
 
+    def get_link_type(self, data_label_1, data_label_2):
+        """Find the type of ``glue`` linking between the given
+        data labels. A link is bi-directional. If there are
+        more than 2 data in the collection, one of the given
+        labels should be the reference data or look-up will fail.
+
+        Parameters
+        ----------
+        data_label_1, data_label_2 : str
+           Labels for the data linked together.
+
+        Returns
+        -------
+        link_type : {'pixels', 'wcs', 'self'}
+            One of the link types accepted by :func:`~jdaviz.configs.imviz.helper.link_image_data`
+            or ``'self'`` if the labels are identical.
+
+        Raises
+        ------
+        ValueError
+            Link look-up failed.
+
+        """
+        if data_label_1 == data_label_2:
+            return "self"
+
+        link_type = None
+        for elink in self.app.data_collection.external_links:
+            elink_labels = (elink.data1.label, elink.data2.label)
+            if data_label_1 in elink_labels and data_label_2 in elink_labels:
+                if isinstance(elink, LinkSame):  # Assumes WCS link never uses LinkSame
+                    link_type = 'pixels'
+                else:  # If not pixels, must be WCS
+                    link_type = 'wcs'
+                break  # Might have duplicate, just grab first match
+
+        if link_type is None:
+            raise ValueError(f'{data_label_1} and {data_label_2} combo not found '
+                             'in data collection external links')
+
+        return link_type
+
     def get_aperture_photometry_results(self):
         """Return aperture photometry results, if any.
         Results are calculated using :ref:`aper-phot-simple` plugin.

jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py

@@ -25,6 +25,8 @@
 
 __all__ = ['SimpleAperturePhotometry']
 
+ASTROPY_LT_5_2 = Version(astropy.__version__) < Version('5.2')
+
 
 @tray_registry('imviz-aper-phot-simple', label="Imviz Simple Aperture Photometry")
 class SimpleAperturePhotometry(PluginTemplateMixin, DatasetSelectMixin, TableMixin, PlotMixin):
@@ -166,9 +168,15 @@ def _subset_selected_changed(self, event={}):
             return
 
         try:
-            self._selected_subset = _get_region_from_spatial_subset(self, subset_selected)
+            self._selected_subset = _get_region_from_spatial_subset(
+                self, self.subset.selected_subset_state)
             self._selected_subset.meta['label'] = subset_selected
-            self.subset_area = int(np.ceil(self._selected_subset.area))
+
+            # Sky subset does not have area. Not worth it to calculate just for a warning.
+            if hasattr(self._selected_subset, 'area'):
+                self.subset_area = int(np.ceil(self._selected_subset.area))
+            else:
+                self.subset_area = 0
 
         except Exception as e:
             self._selected_subset = None
@@ -183,6 +191,8 @@ def _calc_bg_subset_median(self, reg):
         # except here we only care about one stat for the background.
         data = self._selected_data
         comp = data.get_component(data.main_components[0])
+        if hasattr(reg, 'to_pixel'):
+            reg = reg.to_pixel(data.coords)
         aper_mask_stat = reg.to_mask(mode='center')
         img_stat = aper_mask_stat.get_values(comp.data, mask=None)
 
@@ -197,7 +207,7 @@ def _bg_subset_selected_changed(self, event={}):
             return
 
         try:
-            reg = _get_region_from_spatial_subset(self, bg_subset_selected)
+            reg = _get_region_from_spatial_subset(self, self.bg_subset.selected_subset_state)
             self.background_value = self._calc_bg_subset_median(reg)
         except Exception as e:
             self.background_value = 0
@@ -212,8 +222,6 @@ def vue_do_aper_phot(self, *args, **kwargs):
 
         data = self._selected_data
         reg = self._selected_subset
-        xcenter = reg.center.x
-        ycenter = reg.center.y
 
         # Reset last fitted model
         fit_model = None
@@ -226,10 +234,18 @@ def vue_do_aper_phot(self, *args, **kwargs):
                 bg = float(self.background_value)
             except ValueError:  # Clearer error message
                 raise ValueError('Missing or invalid background value')
-            if data.coords is not None:
-                sky_center = data.coords.pixel_to_world(xcenter, ycenter)
+
+            if hasattr(reg, 'to_pixel'):
+                sky_center = reg.center
+                xcenter, ycenter = data.coords.world_to_pixel(sky_center)
             else:
-                sky_center = None
+                xcenter = reg.center.x
+                ycenter = reg.center.y
+                if data.coords is not None:
+                    sky_center = data.coords.pixel_to_world(xcenter, ycenter)
+                else:
+                    sky_center = None
+
             aperture = regions2aperture(reg)
             include_pixarea_fac = False
             include_counts_fac = False
@@ -359,7 +375,7 @@ def vue_do_aper_phot(self, *args, **kwargs):
                     gs = Gaussian1D(amplitude=y_max, mean=x_mean, stddev=std,
                                     fixed={'amplitude': True},
                                     bounds={'amplitude': (y_max * 0.5, y_max)})
-                    if Version(astropy.__version__) < Version('5.2'):
+                    if ASTROPY_LT_5_2:
                         fitter_kw = {}
                     else:
                         fitter_kw = {'filter_non_finite': True}
@@ -523,6 +539,9 @@ def _curve_of_growth(data, centroid, aperture, final_sum, wcs=None, background=0
     """
     n_datapoints += 1  # n + 1
 
+    if hasattr(aperture, 'to_pixel'):
+        aperture = aperture.to_pixel(wcs)
+
     if isinstance(aperture, CircularAperture):
         x_label = 'Radius (pix)'
         x_arr = np.linspace(0, aperture.r, num=n_datapoints)[1:]

jdaviz/configs/imviz/plugins/viewers.py

@@ -3,7 +3,6 @@
 import astropy.units as u
 from astropy.wcs.utils import pixel_to_pixel
 from astropy.visualization import ImageNormalize, LinearStretch, PercentileInterval
-from glue.core.link_helpers import LinkSame
 from glue_jupyter.bqplot.image import BqplotImageView
 
 from jdaviz.configs.imviz import wcs_utils
@@ -281,26 +280,15 @@ def get_link_type(self, data_label):
         if len(self.session.application.data_collection) == 0:
             raise ValueError('No reference data for link look-up')
 
-        # the original links were created against data_collection[0], not necessarily
+        # TODO: Brett Morris might want to look at this for
+        # https://github.com/spacetelescope/jdaviz/pull/2179
+        #     ref_label = self.state.reference_data ???
+        #
+        # The original links were created against data_collection[0], not necessarily
         # against the current viewer reference_data
         ref_label = self.session.application.data_collection[0].label
-        if data_label == ref_label:
-            return 'self'
-
-        link_type = None
-        for elink in self.session.application.data_collection.external_links:
-            elink_labels = (elink.data1.label, elink.data2.label)
-            if data_label in elink_labels and ref_label in elink_labels:
-                if isinstance(elink, LinkSame):  # Assumes WCS link never uses LinkSame
-                    link_type = 'pixels'
-                else:  # If not pixels, must be WCS
-                    link_type = 'wcs'
-                break  # Might have duplicate, just grab first match
-
-        if link_type is None:
-            raise ValueError(f'{data_label} not found in data collection external links')
-
-        return link_type
+
+        return self.jdaviz_helper.get_link_type(ref_label, data_label)
 
     def _get_center_skycoord(self, data=None):
         if data is None: