Fix aperture photometry result for WCS linked dithered non-reference data

CHANGES.rst

@@ -46,6 +46,9 @@ Cubeviz
 Imviz
 ^^^^^
 
+- Fixed inaccurate aperture photometry results when aperture photometry is done on
+  a non-reference image if images are linked by WCS. [#1524]
+
 Mosviz
 ^^^^^^
 

docs/imviz/plugins.rst

@@ -75,6 +75,10 @@ performant at the cost of accuracy but should be accurate to within a pixel
 for most cases. If approximation fails, WCS linking still automatically
 falls back to full transformation.
 
+For the best experience, it is recommended that you decide what kind of
+link you want and set it at the beginning of your Imviz session,
+rather than later.
+
 For more details on linking, see :ref:`dev_glue_linking`.
 
 From the API
@@ -114,12 +118,6 @@ This plugin only considers pixel locations, not sky coordinates.
 Simple Aperture Photometry
 ==========================
 
-.. warning::
-
-    Results for dithered data linked by WCS might be inaccurate unless the selected
-    data is the reference data. See https://github.com/glue-viz/glue-astronomy/issues/52
-    for more details.
-
 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/imviz/helper.py

@@ -500,10 +500,6 @@ def link_image_data(app, link_type='pixels', wcs_fallback_scheme='pixels', wcs_u
     if update_plugin and 'imviz-links-control' in [item['name'] for item in app.state.tray_items]:
         link_plugin = app.get_tray_item_from_name('imviz-links-control')
         link_plugin.linking_in_progress = True
-        app.hub.broadcast(LinkUpdatedMessage(link_type,
-                                             wcs_fallback_scheme == 'pixels',
-                                             wcs_use_affine,
-                                             sender=app))
     else:
         link_plugin = None
 
@@ -573,5 +569,10 @@ def link_image_data(app, link_type='pixels', wcs_fallback_scheme='pixels', wcs_u
             'Images successfully relinked', color='success', timeout=8000, sender=app))
 
     if link_plugin is not None:
+        # Only broadcast after success.
+        app.hub.broadcast(LinkUpdatedMessage(link_type,
+                                             wcs_fallback_scheme == 'pixels',
+                                             wcs_use_affine,
+                                             sender=app))
         # reset the progress spinner
         link_plugin.linking_in_progress = False

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

@@ -2,6 +2,7 @@
 import warnings
 from datetime import datetime
 
+import astropy
 import bqplot
 import numpy as np
 from astropy import units as u
@@ -12,14 +13,15 @@
 from astropy.time import Time
 from glue.core.message import SubsetUpdateMessage
 from ipywidgets import widget_serialization
+from packaging.version import Version
 from photutils.aperture import (ApertureStats, CircularAperture, EllipticalAperture,
                                 RectangularAperture)
 from regions import (CircleAnnulusPixelRegion, CirclePixelRegion, EllipsePixelRegion,
                      RectanglePixelRegion)
 from traitlets import Any, Bool, List, Unicode, observe
 
 from jdaviz.core.custom_traitlets import FloatHandleEmpty
-from jdaviz.core.events import SnackbarMessage
+from jdaviz.core.events import SnackbarMessage, LinkUpdatedMessage
 from jdaviz.core.region_translators import regions2aperture
 from jdaviz.core.registries import tray_registry
 from jdaviz.core.template_mixin import PluginTemplateMixin, DatasetSelectMixin, SubsetSelect
@@ -74,6 +76,7 @@ def __init__(self, *args, **kwargs):
         self._fitted_model_name = 'phot_radial_profile'
 
         self.session.hub.subscribe(self, SubsetUpdateMessage, handler=self._on_subset_update)
+        self.session.hub.subscribe(self, LinkUpdatedMessage, handler=self._on_link_update)
 
     def reset_results(self):
         self.result_available = False
@@ -143,9 +146,20 @@ def _get_region_from_subset(self, subset):
             if subset_grp.label == subset:
                 for sbst in subset_grp.subsets:
                     if sbst.data.label == self.dataset_selected:
-                        # TODO: https://github.com/glue-viz/glue-astronomy/issues/52
-                        return sbst.data.get_selection_definition(
-                                subset_id=subset, format='astropy-regions')
+                        reg = sbst.data.get_selection_definition(
+                            subset_id=subset, format='astropy-regions')
+                        # Works around https://github.com/glue-viz/glue-astronomy/issues/52
+                        # Assume it is always pixel region, not sky region. 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 = self.app._jdaviz_helper.default_viewer
+
+                        x, y, _ = viewer._get_real_xy(
+                            self.app.data_collection[self.dataset_selected],
+                            reg.center.x, reg.center.y)
+                        reg.center.x = x
+                        reg.center.y = y
+                        return reg
         else:
             raise ValueError(f'Subset "{subset}" not found')
 
@@ -159,6 +173,13 @@ def _on_subset_update(self, msg):
         elif sbst.label == self.bg_subset_selected and sbst.data.label == self.dataset_selected:
             self._bg_subset_selected_changed()
 
+    def _on_link_update(self, msg):
+        if self.dataset_selected == '' or self.subset_selected == '':
+            return
+
+        # Force background auto-calculation (including annulus) to update when linking has changed.
+        self._subset_selected_changed()
+
     @observe('subset_selected')
     def _subset_selected_changed(self, event={}):
         subset_selected = event.get('new', self.subset_selected)
@@ -396,8 +417,12 @@ def vue_do_aper_phot(self, *args, **kwargs):
                     gs = Gaussian1D(amplitude=y_max, mean=0, stddev=std,
                                     fixed={'mean': True, 'amplitude': True},
                                     bounds={'amplitude': (y_max * 0.5, y_max)})
+                    if Version(astropy.__version__) <= Version('5.1'):
+                        fitter_kw = {}
+                    else:
+                        fitter_kw = {'filter_non_finite': True}
                     with warnings.catch_warnings(record=True) as warns:
-                        fit_model = fitter(gs, x_data, y_data)
+                        fit_model = fitter(gs, x_data, y_data, **fitter_kw)
                     if len(warns) > 0:
                         msg = os.linesep.join([str(w.message) for w in warns])
                         self.hub.broadcast(SnackbarMessage(
@@ -485,7 +510,13 @@ def _radial_profile(radial_cutout, reg_bb, centroid, raw=False):
     reg_ogrid = np.ogrid[reg_bb.iymin:reg_bb.iymax, reg_bb.ixmin:reg_bb.ixmax]
     radial_dx = reg_ogrid[1] - centroid[0]
     radial_dy = reg_ogrid[0] - centroid[1]
-    radial_r = np.hypot(radial_dx, radial_dy)[~radial_cutout.mask].ravel()  # pix
+    radial_r = np.hypot(radial_dx, radial_dy)
+
+    # Sometimes the mask is smaller than radial_r
+    if radial_cutout.shape != reg_bb.shape:
+        radial_r = radial_r[:radial_cutout.shape[0], :radial_cutout.shape[1]]
+
+    radial_r = radial_r[~radial_cutout.mask].ravel()  # pix
     radial_img = radial_cutout.compressed()  # data unit
 
     if raw:
@@ -495,7 +526,7 @@ def _radial_profile(radial_cutout, reg_bb, centroid, raw=False):
     else:
         # This algorithm is from the imexam package,
         # see licenses/IMEXAM_LICENSE.txt for more details
-        radial_r = list(radial_r)
+        radial_r = np.rint(radial_r).astype(int)
         y_arr = np.bincount(radial_r, radial_img) / np.bincount(radial_r)
         x_arr = np.arange(y_arr.size)
 

jdaviz/configs/imviz/tests/test_simple_aper_phot.py

@@ -73,7 +73,7 @@ def test_plugin_wcs_dithered(self):
             'data_label', 'subset_label', 'timestamp']
         assert_array_equal(tbl['id'], [1, 2])
         assert_allclose(tbl['background'], 0)
-        assert_quantity_allclose(tbl['sum_aper_area'], 63.617251 * (u.pix * u.pix))
+        assert_quantity_allclose(tbl['sum_aper_area'], [63.617251, 62.22684693104279] * (u.pix * u.pix))  # noqa
         assert_array_equal(tbl['pixarea_tot'], None)
         assert_array_equal(tbl['aperture_sum_counts'], None)
         assert_array_equal(tbl['aperture_sum_counts_err'], None)
@@ -99,15 +99,14 @@ def test_plugin_wcs_dithered(self):
         assert_array_equal(tbl['subset_label'], ['Subset 1', 'Subset 1'])
         assert tbl['timestamp'].scale == 'utc'
 
-        # BUG: https://github.com/glue-viz/glue-astronomy/issues/52
-        # Sky should have been the same and the pix different, but not until bug is fixed.
-        # The aperture sum might be different too if mask is off limit in second image.
-        assert_quantity_allclose(tbl['xcentroid'], 4.5 * u.pix)
+        # Sky is the same but xcentroid different due to dithering.
+        # The aperture sum is different too because mask is a little off limit in second image.
+        assert_quantity_allclose(tbl['xcentroid'], [4.5, 5.5] * u.pix)
         assert_quantity_allclose(tbl['ycentroid'], 4.5 * u.pix)
         sky = tbl['sky_centroid']
-        assert_allclose(sky.ra.deg, [337.518943, 337.519241])
-        assert_allclose(sky.dec.deg, [-20.832083, -20.832083])
-        assert_allclose(tbl['sum'], 63.61725123519332)
+        assert_allclose(sky.ra.deg, 337.518943)
+        assert_allclose(sky.dec.deg, -20.832083)
+        assert_allclose(tbl['sum'], [63.617251, 62.22684693104279])
 
         # Make sure it also works on an ellipse subset.
         self.imviz._apply_interactive_region('bqplot:ellipse', (0, 0), (9, 4))