Merge pull request #2240 from larrybradley/aperture-wcs

Improve aperture pixel/sky conversions

Author: larrybradley

CHANGES.rst

@@ -17,6 +17,14 @@ General
 New Features
 ^^^^^^^^^^^^
 
+- ``photutils.aperture``
+
+  - The ``to_sky`` and ``to_pixel`` methods for all pixel and sky
+    aperture types now use the local WCS Jacobian to accurately compute
+    pixel scale factors and rotation angles. This correctly handles WCS
+    with distortions (e.g., SIP polynomial corrections) and non-square
+    pixels. [#2240]
+
 - ``photutils.background``
 
   - Added a ``to_aperture`` method to ``LocalBackground``. [#2118]

docs/user_guide/aperture.rst

@@ -132,6 +132,17 @@ sky apertures, e.g.,::
     >>> pix_aperture  # doctest: +FLOAT_CMP
     <CircularAperture([26.14628817, 56.58410628], r=4.000000000439743)>
 
+.. note::
+
+    Aperture objects require scalar shape parameters (e.g., radius,
+    semi-axes, angle), so only a single reference position can be
+    used for computing local WCS properties (pixel scale, rotation
+    angle). For apertures with multiple positions, the first position
+    is used. Apertures with multiple positions used with a WCS that
+    has spatially-varying distortions may produce inaccurate shape
+    conversions for positions far from the first position.
+
+
 Performing Aperture Photometry
 ------------------------------
 

docs/whats_new/3.0.rst

@@ -414,6 +414,34 @@ parameters::
     ...                         centroid_func=centroid_func)
 
 
+Improved Aperture Pixel-to-Sky Conversions
+==========================================
+
+The :meth:`~photutils.aperture.PixelAperture.to_sky` and
+:meth:`~photutils.aperture.SkyAperture.to_pixel` methods for all
+aperture types have been significantly improved. They now use the local
+WCS Jacobian to accurately compute pixel scale factors and rotation
+angles, correctly handling WCS with distortions (e.g., SIP polynomial
+corrections) and non-square pixels.
+
+For elliptical and rectangular apertures specifically, the conversion
+now uses singular value decomposition (SVD) of the local Jacobian to
+compute independent scale factors along the width and height axes. This
+is more accurate than applying a single mean pixel scale to all linear
+dimensions, particularly near the edge of a distorted WCS or when the
+pixel scale differs significantly along the two axes.
+
+For all aperture types, the WCS properties (pixel scale, rotation angle)
+are evaluated at the first aperture position. Because aperture objects
+require scalar shape parameters, only a single reference position can
+be used. For apertures with multiple positions used with a WCS that
+has spatially-varying distortions, this may produce inaccurate shape
+conversions for positions far from the first position. This behavior
+is now clearly documented in the :ref:`photutils-aperture` user guide
+and in the ``Notes`` section of each ``to_sky`` and ``to_pixel`` method
+docstring.
+
+
 API Changes
 ===========
 

photutils/aperture/circle.py

@@ -6,6 +6,9 @@
 
 import math
 
+import astropy.units as u
+import numpy as np
+from astropy.coordinates import Angle
 from astropy.utils import lazyproperty
 
 from photutils.aperture.attributes import (PixelPositions, PositiveScalar,
@@ -15,6 +18,8 @@
 from photutils.aperture.mask import ApertureMask
 from photutils.geometry import circular_overlap_grid
 from photutils.utils._deprecation import deprecated_positional_kwargs
+from photutils.utils._wcs_helpers import (pixel_to_sky_mean_scale,
+                                          sky_to_pixel_mean_scale)
 
 __all__ = [
     'CircularAnnulus',
@@ -222,8 +227,26 @@ def to_sky(self, wcs):
         -------
         aperture : `SkyCircularAperture` object
             A `SkyCircularAperture` object.
+
+        Notes
+        -----
+        The aperture shape parameters are converted using the local
+        WCS pixel scale evaluated at the first aperture position.
+        Because aperture objects require scalar shape parameters, only
+        a single reference position is used for the conversion. For
+        apertures with multiple positions used with a WCS that has
+        spatially-varying distortions, this may produce inaccurate
+        results for positions far from the first position.
         """
-        return SkyCircularAperture(**self._to_sky_params(wcs))
+        xpos, ypos = np.transpose(self.positions)
+        positions = wcs.pixel_to_world(xpos, ypos)
+
+        first_pos = np.atleast_2d(self.positions)[0]
+        _, mean_scale = pixel_to_sky_mean_scale(
+            (float(first_pos[0]), float(first_pos[1])), wcs)
+
+        r = Angle(self.r * mean_scale, 'arcsec')
+        return SkyCircularAperture(positions=positions, r=r)
 
 
 class CircularAnnulus(CircularMaskMixin, PixelAperture):
@@ -356,8 +379,27 @@ def to_sky(self, wcs):
         -------
         aperture : `SkyCircularAnnulus` object
             A `SkyCircularAnnulus` object.
+
+        Notes
+        -----
+        The aperture shape parameters are converted using the local
+        WCS pixel scale evaluated at the first aperture position.
+        Because aperture objects require scalar shape parameters, only
+        a single reference position is used for the conversion. For
+        apertures with multiple positions used with a WCS that has
+        spatially-varying distortions, this may produce inaccurate
+        results for positions far from the first position.
         """
-        return SkyCircularAnnulus(**self._to_sky_params(wcs))
+        xpos, ypos = np.transpose(self.positions)
+        positions = wcs.pixel_to_world(xpos, ypos)
+
+        first_pos = np.atleast_2d(self.positions)[0]
+        _, mean_scale = pixel_to_sky_mean_scale(
+            (float(first_pos[0]), float(first_pos[1])), wcs)
+
+        r_in = Angle(self.r_in * mean_scale, 'arcsec')
+        r_out = Angle(self.r_out * mean_scale, 'arcsec')
+        return SkyCircularAnnulus(positions=positions, r_in=r_in, r_out=r_out)
 
 
 class SkyCircularAperture(SkyAperture):
@@ -410,8 +452,25 @@ def to_pixel(self, wcs):
         -------
         aperture : `CircularAperture` object
             A `CircularAperture` object.
+
+        Notes
+        -----
+        The aperture shape parameters are converted using the local
+        WCS pixel scale evaluated at the first aperture position.
+        Because aperture objects require scalar shape parameters, only
+        a single reference position is used for the conversion. For
+        apertures with multiple positions used with a WCS that has
+        spatially-varying distortions, this may produce inaccurate
+        results for positions far from the first position.
         """
-        return CircularAperture(**self._to_pixel_params(wcs))
+        xpos, ypos = wcs.world_to_pixel(self.positions)
+        positions = np.transpose((xpos, ypos))
+
+        skypos = self.positions if self.isscalar else self.positions[0]
+        _, mean_scale = sky_to_pixel_mean_scale(skypos, wcs)
+
+        r = self.r.to(u.arcsec).value * mean_scale
+        return CircularAperture(positions=positions, r=r)
 
 
 class SkyCircularAnnulus(SkyAperture):
@@ -473,5 +532,23 @@ def to_pixel(self, wcs):
         -------
         aperture : `CircularAnnulus` object
             A `CircularAnnulus` object.
+
+        Notes
+        -----
+        The aperture shape parameters are converted using the local
+        WCS pixel scale evaluated at the first aperture position.
+        Because aperture objects require scalar shape parameters, only
+        a single reference position is used for the conversion. For
+        apertures with multiple positions used with a WCS that has
+        spatially-varying distortions, this may produce inaccurate
+        results for positions far from the first position.
         """
-        return CircularAnnulus(**self._to_pixel_params(wcs))
+        xpos, ypos = wcs.world_to_pixel(self.positions)
+        positions = np.transpose((xpos, ypos))
+
+        skypos = self.positions if self.isscalar else self.positions[0]
+        _, mean_scale = sky_to_pixel_mean_scale(skypos, wcs)
+
+        r_in = self.r_in.to(u.arcsec).value * mean_scale
+        r_out = self.r_out.to(u.arcsec).value * mean_scale
+        return CircularAnnulus(positions=positions, r_in=r_in, r_out=r_out)

photutils/aperture/core.py

@@ -8,14 +8,12 @@
 import warnings
 from copy import deepcopy
 
-import astropy.units as u
 import numpy as np
 from astropy.coordinates import SkyCoord
 from astropy.utils import lazyproperty
 
 from photutils.aperture.bounding_box import BoundingBox
 from photutils.utils._deprecation import deprecated_positional_kwargs
-from photutils.utils._wcs_helpers import _pixel_scale_angle_at_skycoord
 
 __all__ = ['Aperture', 'PixelAperture', 'SkyAperture']
 
@@ -714,59 +712,6 @@ def plot(self, ax=None, origin=(0, 0), **kwargs):
 
         return patches
 
-    def _to_sky_params(self, wcs):
-        """
-        Convert the pixel aperture parameters to those for a sky
-        aperture.
-
-        Parameters
-        ----------
-        wcs : WCS object
-            A world coordinate system (WCS) transformation that
-            supports the `astropy shared interface for WCS
-            <https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_
-            (e.g., `astropy.wcs.WCS`, `gwcs.wcs.WCS`).
-
-        Returns
-        -------
-        sky_params : dict
-            A dictionary of parameters for an equivalent sky aperture.
-        """
-        sky_params = {}
-        xpos, ypos = np.transpose(self.positions)
-        sky_params['positions'] = skypos = wcs.pixel_to_world(xpos, ypos)
-
-        # Aperture objects require scalar shape parameters (e.g.,
-        # radius, a, b, theta, etc.), therefore we must calculate the
-        # pixel scale and angle at only a single sky position, which
-        # we take as the first aperture position. For apertures with
-        # multiple positions used with a WCS that contains distortions
-        # (e.g., a spatially-dependent pixel scale), this may lead to
-        # unexpected results (e.g., results that are dependent of the
-        # order of the positions). There is no good way to fix this with
-        # the current Aperture API allowing multiple positions.
-        if not self.isscalar:
-            skypos = skypos[0]
-        _, pixscale, angle = _pixel_scale_angle_at_skycoord(skypos, wcs)
-
-        for param in self._params:
-            value = getattr(self, param)
-            if param == 'positions':
-                continue
-
-            if param == 'theta':
-                # photutils aperture sky angles are defined as the PA of
-                # the semimajor axis (i.e., relative to the WCS latitude
-                # axis). region sky angles are defined relative to the WCS
-                # longitude axis.
-                value = value - angle.to(u.rad)
-            else:
-                value = (value * u.pix * pixscale).to(u.arcsec)
-
-            sky_params[param] = value
-
-        return sky_params
-
     @abc.abstractmethod
     def to_sky(self, wcs):
         """
@@ -794,61 +739,6 @@ class SkyAperture(Aperture):
     coordinates.
     """
 
-    def _to_pixel_params(self, wcs):
-        """
-        Convert the sky aperture parameters to those for a pixel
-        aperture.
-
-        Parameters
-        ----------
-        wcs : WCS object
-            A world coordinate system (WCS) transformation that
-            supports the `astropy shared interface for WCS
-            <https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_
-            (e.g., `astropy.wcs.WCS`, `gwcs.wcs.WCS`).
-
-        Returns
-        -------
-        pixel_params : dict
-            A dictionary of parameters for an equivalent pixel aperture.
-        """
-        pixel_params = {}
-
-        xpos, ypos = wcs.world_to_pixel(self.positions)
-        pixel_params['positions'] = np.transpose((xpos, ypos))
-
-        # Aperture objects require scalar shape parameters (e.g.,
-        # radius, a, b, theta, etc.), therefore we must calculate the
-        # pixel scale and angle at only a single sky position, which
-        # we take as the first aperture position. For apertures with
-        # multiple positions used with a WCS that contains distortions
-        # (e.g., a spatially-dependent pixel scale), this may lead to
-        # unexpected results (e.g., results that are dependent of the
-        # order of the positions). There is no good way to fix this with
-        # the current Aperture API allowing multiple positions.
-        skypos = self.positions if self.isscalar else self.positions[0]
-        _, pixscale, angle = _pixel_scale_angle_at_skycoord(skypos, wcs)
-
-        for param in self._params:
-            value = getattr(self, param)
-            if param == 'positions':
-                continue
-
-            if param == 'theta':
-                # photutils aperture sky angles are defined as the PA of
-                # the semimajor axis (i.e., relative to the WCS latitude
-                # axis). region sky angles are defined relative to the WCS
-                # longitude axis.
-                value = (value + angle).to(u.radian)
-            elif value.unit.physical_type == 'angle':
-                value = (value / pixscale).to(u.pixel).value
-            else:
-                value = value.value
-
-            pixel_params[param] = value
-
-        return pixel_params
-
     @abc.abstractmethod
     def to_pixel(self, wcs):
         """