Moar fixes and tests

CHANGES.rst

@@ -41,6 +41,9 @@ Imviz
   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
 ^^^^^^
 

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

@@ -357,29 +357,39 @@ 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, wcs=data.coords)
 
-            # 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)
+            # 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:
             self.set_center(self.get_center(), update=False)
             self.hub.broadcast(SnackbarMessage(
-                f"Failed to calculate centroid: {repr(err)}", color='error', sender=self))
+                f"Failed to calculate centroid: {err!r}", color='error', sender=self))
         else:
             self.set_center((x, y), update=True)
 

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
+        label 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

@@ -168,7 +168,8 @@ 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
 
             # Sky subset does not have area. Not worth it to calculate just for a warning.
@@ -206,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
@@ -246,6 +247,13 @@ def vue_do_aper_phot(self, *args, **kwargs):
                     sky_center = None
 
             aperture = regions2aperture(reg)
+
+# UNTIL HERE - pass in pixel region explicitly and see if answers are correct
+#            print("HERE 1", reg)
+#            print("HERE 2", aperture)
+#            if hasattr(reg, "to_pixel"):
+#                print("HERE 3", reg.to_pixel(data.coords))
+
             include_pixarea_fac = False
             include_counts_fac = False
             include_flux_scale = False

jdaviz/configs/imviz/plugins/viewers.py

@@ -2,7 +2,6 @@
 
 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
@@ -280,23 +279,12 @@ 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)

jdaviz/configs/imviz/tests/test_simple_aper_phot.py

@@ -7,7 +7,8 @@
 from numpy.testing import assert_allclose, assert_array_equal
 from photutils.aperture import (ApertureStats, CircularAperture, EllipticalAperture,
                                 RectangularAperture, EllipticalAnnulus)
-from regions import CircleAnnulusPixelRegion, CirclePixelRegion, RectanglePixelRegion, PixCoord
+from regions import (CircleAnnulusPixelRegion, CirclePixelRegion, EllipsePixelRegion,
+                     RectanglePixelRegion, PixCoord)
 
 from jdaviz.configs.imviz.plugins.aper_phot_simple.aper_phot_simple import (
     _curve_of_growth, _radial_profile)
@@ -206,10 +207,13 @@ def setup_class(self, imviz_helper):
 
         # Regions to be used for aperture photometry
         regions = []
-        positions = [(145.1, 168.3), (84.7, 224.1), (48.3, 200.3)]
+        positions = [(145.1, 168.3), (48.3, 200.3)]
         for x, y in positions:
             regions.append(CirclePixelRegion(center=PixCoord(x=x, y=y), radius=5))
-        regions.append(RectanglePixelRegion(center=PixCoord(x=229, y=152), width=17, height=5))
+        regions += [
+            EllipsePixelRegion(center=PixCoord(x=84.7, y=224.1), width=23, height=9,
+                               angle=2.356 * u.rad),
+            RectanglePixelRegion(center=PixCoord(x=229, y=152), width=17, height=7)]
         imviz_helper.load_regions(regions)
 
         self.imviz = imviz_helper
@@ -224,9 +228,9 @@ def setup_class(self, imviz_helper):
         ('no_wcs', 5.0)])
     @pytest.mark.parametrize(('subset_label', 'expected_sum'), [
         ('Subset 1', 738.8803424408962),
-        ('Subset 2', 348.33262363800094),
-        ('Subset 3', 857.5194857987592),
-        ('Subset 4', 762.3263959884644)])
+        ('Subset 2', 857.5194857987592),
+        ('Subset 3', 472.17364321556005),
+        ('Subset 4', 837.0023608207703)])
     def test_aperphot(self, data_label, local_bkg, subset_label, expected_sum):
         """All data should give similar result for the same Subset."""
         self.phot_plugin.dataset_selected = data_label
@@ -235,7 +239,6 @@ def test_aperphot(self, data_label, local_bkg, subset_label, expected_sum):
         self.phot_plugin.vue_do_aper_phot()
         tbl = self.imviz.get_aperture_photometry_results()
 
-        # FIXME: Subset 4 (rectangle) has huge difference, why???
         # Imperfect down-sampling and imperfect apertures, so 10% is good enough.
         assert_allclose(tbl['sum'][-1], expected_sum, rtol=0.1)
 

jdaviz/core/region_translators.py

@@ -20,8 +20,8 @@
 __all__ = ['regions2roi', 'regions2aperture', 'aperture2regions']
 
 
-def _get_region_from_spatial_subset(plugin_obj, subset_label):
-    """Convert a ``glue`` ROI in the given subset to ``regions`` shape.
+def _get_region_from_spatial_subset(plugin_obj, subset_state):
+    """Convert the given ``glue`` ROI subset state to ``regions`` shape.
 
     .. note:: This is for internal use only in Imviz plugins.
 
@@ -31,50 +31,37 @@ def _get_region_from_spatial_subset(plugin_obj, subset_label):
         Plugin instance that needs this translation.
         The plugin is assumed to have a special setup that gives
         it access to these attributes: ``app`` and ``dataset_selected``.
-        The ``app._jdaviz_helper.default_viewer`` attribute must also
-        exist and point to an image viewer that has a ``get_link_type``
-        method.
+        The ``app._jdaviz_helper.get_link_type`` method must also
+        exist.
 
-    subset_label : str
-        Label of selected subset in the plugin.
+    subset_state : obj
+        ROI subset state to translate.
 
     Returns
     -------
     reg : `regions.Region`
         An equivalent ``regions`` shape. This can be a pixel or sky
         region, so the plugin needs to be able to deal with both.
 
-    Raises
-    ------
-    ValueError
-        Subset is not found.
-
     See Also
     --------
     regions2roi
 
     """
     from glue_astronomy.translators.regions import roi_subset_state_to_region
 
-    for subset_grp in plugin_obj.app.data_collection.subset_groups:
-        if subset_grp.label == subset_label:
-            sbst = subset_grp.subset_state
+    # Subset is defined against its parent. This is not necessarily
+    # the current viewer reference data for when we can change that.
+    # See https://github.com/spacetelescope/jdaviz/issues/2230
+    link_type = plugin_obj.app._jdaviz_helper.get_link_type(
+        subset_state.xatt.parent.label, plugin_obj.dataset_selected)
 
-            # Works around https://github.com/glue-viz/glue-astronomy/issues/52
-            # Even with multiple viewers, they all seem to share the same
-            # reference image even when it is not loaded in all the viewers,
-            # so use default viewer.
-            viewer = plugin_obj.app._jdaviz_helper.default_viewer
-            link_type = viewer.get_link_type(plugin_obj.dataset_selected)
-
-            if link_type == 'wcs':
-                to_sky = True
-            else:
-                to_sky = False
-
-            return roi_subset_state_to_region(sbst, to_sky=to_sky)
+    if link_type == 'wcs':
+        to_sky = True
+    else:
+        to_sky = False
 
-    raise ValueError(f'Subset "{subset_label}" not found')
+    return roi_subset_state_to_region(subset_state, to_sky=to_sky)
 
 
 def regions2roi(region_shape, wcs=None):
@@ -200,7 +187,7 @@ def regions2aperture(region_shape):
     elif isinstance(region_shape, EllipseSkyRegion):
         aperture = SkyEllipticalAperture(
             region_shape.center, region_shape.width * 0.5, region_shape.height * 0.5,
-            theta=region_shape.angle)
+            theta=(region_shape.angle - (90 * u.deg)))
 
     elif isinstance(region_shape, RectanglePixelRegion):
         aperture = RectangularAperture(
@@ -209,7 +196,8 @@ def regions2aperture(region_shape):
 
     elif isinstance(region_shape, RectangleSkyRegion):
         aperture = SkyRectangularAperture(
-            region_shape.center, region_shape.width, region_shape.height, theta=region_shape.angle)
+            region_shape.center, region_shape.width, region_shape.height,
+            theta=(region_shape.angle - (90 * u.deg)))
 
     elif isinstance(region_shape, CircleAnnulusPixelRegion):
         aperture = CircularAnnulus(
@@ -229,7 +217,7 @@ def regions2aperture(region_shape):
         aperture = SkyEllipticalAnnulus(
             region_shape.center, region_shape.inner_width * 0.5, region_shape.outer_width * 0.5,
             region_shape.outer_height * 0.5, b_in=region_shape.inner_height * 0.5,
-            theta=region_shape.angle)
+            theta=(region_shape.angle - (90 * u.deg)))
 
     elif isinstance(region_shape, RectangleAnnulusPixelRegion):
         aperture = RectangularAnnulus(
@@ -240,7 +228,8 @@ def regions2aperture(region_shape):
     elif isinstance(region_shape, RectangleAnnulusSkyRegion):
         aperture = SkyRectangularAnnulus(
             region_shape.center, region_shape.inner_width, region_shape.outer_width,
-            region_shape.outer_height, h_in=region_shape.inner_height, theta=region_shape.angle)
+            region_shape.outer_height, h_in=region_shape.inner_height,
+            theta=(region_shape.angle - (90 * u.deg)))
 
     else:
         raise NotImplementedError(f'{region_shape.__class__.__name__} is not supported')
@@ -299,7 +288,7 @@ def aperture2regions(aperture):
     elif isinstance(aperture, SkyEllipticalAperture):
         region_shape = EllipseSkyRegion(
             center=aperture.positions, width=aperture.a * 2, height=aperture.b * 2,
-            angle=aperture.theta)
+            angle=(aperture.theta + (90 * u.deg)))
 
     elif isinstance(aperture, RectangularAperture):
         region_shape = RectanglePixelRegion(
@@ -308,7 +297,8 @@ def aperture2regions(aperture):
 
     elif isinstance(aperture, SkyRectangularAperture):
         region_shape = RectangleSkyRegion(
-            center=aperture.positions, width=aperture.w, height=aperture.h, angle=aperture.theta)
+            center=aperture.positions, width=aperture.w, height=aperture.h,
+            angle=(aperture.theta + (90 * u.deg)))
 
     elif isinstance(aperture, CircularAnnulus):
         region_shape = CircleAnnulusPixelRegion(
@@ -329,7 +319,8 @@ def aperture2regions(aperture):
         region_shape = EllipseAnnulusSkyRegion(
             center=aperture.positions, inner_width=aperture.a_in * 2,
             inner_height=aperture.b_in * 2, outer_width=aperture.a_out * 2,
-            outer_height=aperture.b_out * 2, angle=aperture.theta)
+            outer_height=aperture.b_out * 2,
+            angle=(aperture.theta + (90 * u.deg)))
 
     elif isinstance(aperture, RectangularAnnulus):
         region_shape = RectangleAnnulusPixelRegion(
@@ -340,7 +331,8 @@ def aperture2regions(aperture):
     elif isinstance(aperture, SkyRectangularAnnulus):
         region_shape = RectangleAnnulusSkyRegion(
             center=aperture.positions, inner_width=aperture.w_in, inner_height=aperture.h_in,
-            outer_width=aperture.w_out, outer_height=aperture.h_out, angle=aperture.theta)
+            outer_width=aperture.w_out, outer_height=aperture.h_out,
+            angle=(aperture.theta + (90 * u.deg)))
 
     else:  # pragma: no cover
         raise NotImplementedError(f'{aperture.__class__.__name__} is not supported')
@@ -362,5 +354,5 @@ def positions2pixcoord(positions):
 
 
 def theta2angle(theta):
-    """Convert ``photutils`` theta to ``regions`` angle."""
+    """Convert ``photutils`` theta to ``regions`` angle for pixel regions."""
     return Angle(theta, u.radian)