Add use_ellipse keyword to circle to_sky and to_pixel methods

CHANGES.rst

@@ -26,6 +26,14 @@ New Features
   more accurate transformations for directed regions such as ellipses and
   rectangles. [#650, #651]
 
+- Added a ``use_ellipse`` keyword to the ``to_pixel`` and ``to_sky``
+  methods of circular and circular annulus region classes that allows
+  the user to specify whether to return an elliptical region instead of
+  a circular region when converting between pixel and sky coordinates.
+  This is useful for cases where the WCS is distorted or has different
+  pixel scales along different axes and the circular region would be
+  significantly distorted in pixel or sky coordinates. [#652]
+
 - Improved the string representations of all ``Regions`` and ``Region``
   objects. [#653]
 

regions/_utils/wcs_helpers.py

@@ -797,6 +797,136 @@ def sky_ellipse_to_pixel_svd(skycoord, wcs, width_arcsec, height_arcsec,
     return center, pixel_width, pixel_height, pixel_angle
 
 
+def sky_to_pixel_svd_scales(skycoord, wcs):
+    """
+    Compute the pixel center, principal-axis scale factors, and pixel
+    angle for a sky-to-pixel conversion using SVD of the local Jacobian.
+
+    Uses the singular value decomposition (SVD) of the local Jacobian
+    ``J = d(pixel)/d(sky_arcsec)`` to find the natural principal axes
+    of the WCS transformation at the given sky position. The singular
+    values give the scale factors along the major and minor axes of the
+    ellipse that a unit circle in sky space maps to in pixel space. The
+    left singular vectors give the directions of those axes in pixel
+    coordinates.
+
+    This is the appropriate method for converting a circular sky region
+    to a pixel ellipse, as the resulting ellipse accurately represents
+    the true shape of the WCS mapping (i.e., the tightest-fitting pixel
+    ellipse that contains the sky circle).
+
+    Parameters
+    ----------
+    skycoord : `~astropy.coordinates.SkyCoord`
+        The sky coordinate of the region center.
+
+    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
+    -------
+    center : `~regions.PixCoord`
+        The pixel center position.
+
+    scale_major : float
+        The scale factor along the major (maximum-stretch) axis
+        (pixels per arcsec).
+
+    scale_minor : float
+        The scale factor along the minor (minimum-stretch) axis
+        (pixels per arcsec).
+
+    pixel_angle : `~astropy.coordinates.Angle`
+        The pixel rotation angle of the major axis, measured
+        counterclockwise from the positive x-axis, wrapped to
+        [0, 360) degrees.
+    """
+    x0, y0 = wcs.world_to_pixel(skycoord)
+    center = PixCoord(x=x0, y=y0)
+
+    jacobian = compute_local_wcs_jacobian(skycoord, wcs)
+    u_mat, s_vals, _vt = np.linalg.svd(jacobian)
+
+    # Pixel angle of the major axis: direction of u_mat[:,0] in pixel
+    # space. No parity correction needed — pixel space has no axis
+    # reflection.
+    pixel_angle = Angle(
+        np.rad2deg(np.arctan2(u_mat[1, 0], u_mat[0, 0])) * u.deg).wrap_at(
+            360 * u.deg)
+
+    return center, s_vals[0], s_vals[1], pixel_angle
+
+
+def pixel_to_sky_svd_scales(pixcoord, wcs):
+    """
+    Compute the sky center, principal-axis scale factors, and sky angle
+    for a pixel-to-sky conversion using SVD of the inverse Jacobian.
+
+    Uses the singular value decomposition (SVD) of the local inverse
+    Jacobian ``J^{-1} = d(sky)/d(pixel)`` to find the natural principal
+    axes of the WCS transformation at the given pixel position. The
+    singular values give the scale factors along the major and minor
+    axes of the ellipse that a unit circle in pixel space maps to in sky
+    space. The left singular vectors give the directions of those axes
+    in tangent-plane coordinates.
+
+    This is the appropriate method for converting a circular pixel
+    region to a sky ellipse, as the resulting ellipse accurately
+    represents the true shape of the WCS mapping (i.e., the
+    tightest-fitting sky ellipse that contains the pixel circle).
+
+    Parameters
+    ----------
+    pixcoord : `~regions.PixCoord`
+        The pixel coordinate of the region center.
+
+    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
+    -------
+    center : `~astropy.coordinates.SkyCoord`
+        The sky center position.
+
+    scale_major : float
+        The scale factor along the major (maximum-stretch) axis
+        (arcsec per pixel).
+
+    scale_minor : float
+        The scale factor along the minor (minimum-stretch) axis
+        (arcsec per pixel).
+
+    sky_angle : `~astropy.coordinates.Angle`
+        The sky rotation angle of the major axis, measured
+        counterclockwise from the longitude (RA) axis, wrapped to
+        [0, 360) degrees.
+    """
+    center = wcs.pixel_to_world(pixcoord.x, pixcoord.y)
+
+    jacobian = compute_local_wcs_jacobian(center, wcs)
+    jacobian_inv = np.linalg.inv(jacobian)
+
+    # WCS parity: standard WCS has RA increasing to the left, giving
+    # det(J) < 0. Apply parity to the RA component of U[:,0] to recover
+    # the correct sky angle from the tangent-plane direction.
+    parity = np.sign(np.linalg.det(jacobian))
+
+    u_mat, s_vals, _vt = np.linalg.svd(jacobian_inv)
+
+    # Sky angle of the major axis in tangent-plane coordinates
+    sky_angle = Angle(
+        np.rad2deg(np.arctan2(u_mat[1, 0], parity * u_mat[0, 0])) * u.deg
+    ).wrap_at(360 * u.deg)
+
+    return center, s_vals[0], s_vals[1], sky_angle
+
+
 def wcs_pixel_scale_angle(skycoord, wcs):
     """
     Calculate the pixel coordinate, scale, and WCS rotation angle at the
@@ -846,10 +976,10 @@ def wcs_pixel_scale_angle(skycoord, wcs):
     cdelt_y = sky0.separation(sky_y).arcsec
     scale = np.sqrt(cdelt_x * cdelt_y)
 
-    # Compute the angle by offsetting in latitude by exactly the
-    # local cdelt (geometric-mean pixel scale in degrees). This
-    # ensures the finite-difference derivative probes the same scale
-    # of the distortion field.
+    # Compute the angle by offsetting in latitude by exactly the local
+    # cdelt (geometric-mean pixel scale in degrees). This ensures
+    # the finite-difference derivative probes the same scale of the
+    # distortion field.
     cdelt_deg = scale / 3600  # arcsec -> deg
     skycoord_offset = skycoord.directional_offset_by(
         0.0, cdelt_deg * u.deg)

regions/shapes/annulus.py

@@ -11,9 +11,11 @@
 from regions._utils.wcs_helpers import (pixel_ellipse_to_sky_svd,
                                         pixel_to_sky_mean_scale,
                                         pixel_to_sky_scales,
+                                        pixel_to_sky_svd_scales,
                                         sky_ellipse_to_pixel_svd,
                                         sky_to_pixel_mean_scale,
-                                        sky_to_pixel_scales)
+                                        sky_to_pixel_scales,
+                                        sky_to_pixel_svd_scales)
 from regions.core.attributes import (PositiveScalar, PositiveScalarAngle,
                                      RegionMetaDescr, RegionVisualDescr,
                                      ScalarAngle, ScalarPixCoord,
@@ -153,7 +155,43 @@ def _outer_region(self):
         return self._component_class(self.center, self.outer_radius,
                                      self.meta, self.visual)
 
-    def to_sky(self, wcs):
+    def to_sky(self, wcs, *, use_ellipse=False):
+        """
+        Return a sky region from this pixel region.
+
+        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`).
+
+        use_ellipse : bool, optional
+            If `True`, return an `~regions.EllipseAnnulusSkyRegion`
+            instead of a `~regions.CircleAnnulusSkyRegion`. An ellipse
+            annulus is generally a better approximation when the WCS has
+            distortions or different pixel scales along different axes.
+            Default is `False`.
+
+        Returns
+        -------
+        region : `~regions.CircleAnnulusSkyRegion` or `~regions.EllipseAnnulusSkyRegion`
+            The sky region. An ellipse annulus is returned if
+            ``use_ellipse`` is `True`.
+        """
+        if use_ellipse:
+            center, scale_major, scale_minor, angle = pixel_to_sky_svd_scales(
+                self.center, wcs)
+            inner_width = 2 * self.inner_radius * scale_major * u.arcsec
+            outer_width = 2 * self.outer_radius * scale_major * u.arcsec
+            inner_height = 2 * self.inner_radius * scale_minor * u.arcsec
+            outer_height = 2 * self.outer_radius * scale_minor * u.arcsec
+            return EllipseAnnulusSkyRegion(
+                center, inner_width, outer_width,
+                inner_height, outer_height, angle=angle,
+                meta=self.meta.copy(), visual=self.visual.copy())
+
         center, mean_scale = pixel_to_sky_mean_scale(self.center, wcs)
         inner_radius = self.inner_radius * mean_scale * u.arcsec
         outer_radius = self.outer_radius * mean_scale * u.arcsec
@@ -224,7 +262,45 @@ def __init__(self, center, inner_radius, outer_radius, meta=None,
         if inner_radius >= outer_radius:
             raise ValueError('outer_radius must be greater than inner_radius')
 
-    def to_pixel(self, wcs):
+    def to_pixel(self, wcs, *, use_ellipse=False):
+        """
+        Return a pixel region from this sky region.
+
+        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`).
+
+        use_ellipse : bool, optional
+            If `True`, return an `~regions.EllipseAnnulusPixelRegion`
+            instead of a `~regions.CircleAnnulusPixelRegion`. An ellipse
+            annulus is generally a better approximation when the WCS has
+            distortions or different pixel scales along different axes.
+            Default is `False`.
+
+        Returns
+        -------
+        region : `~regions.CircleAnnulusPixelRegion` or `~regions.EllipseAnnulusPixelRegion`
+            The pixel region. An ellipse annulus is returned if
+            ``use_ellipse`` is `True`.
+        """
+        if use_ellipse:
+            center, scale_major, scale_minor, angle = sky_to_pixel_svd_scales(
+                self.center, wcs)
+            inner_radius_arcsec = self.inner_radius.to(u.arcsec).value
+            outer_radius_arcsec = self.outer_radius.to(u.arcsec).value
+            inner_width = 2 * inner_radius_arcsec * scale_major
+            outer_width = 2 * outer_radius_arcsec * scale_major
+            inner_height = 2 * inner_radius_arcsec * scale_minor
+            outer_height = 2 * outer_radius_arcsec * scale_minor
+            return EllipseAnnulusPixelRegion(
+                center, inner_width, outer_width,
+                inner_height, outer_height, angle=angle,
+                meta=self.meta.copy(), visual=self.visual.copy())
+
         center, mean_scale = sky_to_pixel_mean_scale(self.center, wcs)
         inner_radius = self.inner_radius.to(u.arcsec).value * mean_scale
         outer_radius = self.outer_radius.to(u.arcsec).value * mean_scale

regions/shapes/circle.py

@@ -11,7 +11,9 @@
 
 from regions._geometry import circular_overlap_grid
 from regions._utils.wcs_helpers import (pixel_to_sky_mean_scale,
-                                        sky_to_pixel_mean_scale)
+                                        pixel_to_sky_svd_scales,
+                                        sky_to_pixel_mean_scale,
+                                        sky_to_pixel_svd_scales)
 from regions.core.attributes import (PositiveScalar, PositiveScalarAngle,
                                      RegionMetaDescr, RegionVisualDescr,
                                      ScalarPixCoord, ScalarSkyCoord)
@@ -20,6 +22,7 @@
 from regions.core.mask import RegionMask
 from regions.core.metadata import RegionMeta, RegionVisual
 from regions.core.pixcoord import PixCoord
+from regions.shapes.ellipse import EllipsePixelRegion, EllipseSkyRegion
 from regions.shapes.polygon import PolygonPixelRegion
 
 __all__ = ['CirclePixelRegion', 'CircleSkyRegion']
@@ -86,7 +89,40 @@ def contains(self, pixcoord):
         else:
             return np.logical_not(in_circle)
 
-    def to_sky(self, wcs):
+    def to_sky(self, wcs, *, use_ellipse=False):
+        """
+        Return a sky region from this pixel region.
+
+        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`).
+
+        use_ellipse : bool, optional
+            If `True`, return an `~regions.EllipseSkyRegion` instead
+            of a `~regions.CircleSkyRegion`. An ellipse is generally
+            a better approximation when the WCS has distortions or
+            different pixel scales along different axes. Default is
+            `False`.
+
+        Returns
+        -------
+        region : `~regions.CircleSkyRegion` or `~regions.EllipseSkyRegion`
+            The sky region. An ellipse is returned if ``use_ellipse``
+            is `True`.
+        """
+        if use_ellipse:
+            center, scale_major, scale_minor, angle = pixel_to_sky_svd_scales(
+                self.center, wcs)
+            width = Angle(2 * self.radius * scale_major, 'arcsec')
+            height = Angle(2 * self.radius * scale_minor, 'arcsec')
+            return EllipseSkyRegion(center, width, height, angle=angle,
+                                    meta=self.meta.copy(),
+                                    visual=self.visual.copy())
+
         center, mean_scale = pixel_to_sky_mean_scale(self.center, wcs)
         radius = Angle(self.radius * mean_scale, 'arcsec')
         return CircleSkyRegion(center, radius, meta=self.meta.copy(),
@@ -233,7 +269,42 @@ def __init__(self, center, radius, meta=None, visual=None):
         self.meta = meta or RegionMeta()
         self.visual = visual or RegionVisual()
 
-    def to_pixel(self, wcs):
+    def to_pixel(self, wcs, *, use_ellipse=False):
+        """
+        Return a pixel region from this sky region.
+
+        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`).
+
+        use_ellipse : bool, optional
+            If `True`, return an `~regions.EllipsePixelRegion` instead
+            of a `~regions.CirclePixelRegion`. An ellipse is generally
+            a better approximation when the WCS has distortions or
+            different pixel scales along different axes. Default is
+            `False`.
+
+        Returns
+        -------
+        region : `~regions.CirclePixelRegion` or `~regions.EllipsePixelRegion`
+            The pixel region. An ellipse is returned if ``use_ellipse``
+            is `True`.
+        """
+        if use_ellipse:
+            center, scale_major, scale_minor, angle = sky_to_pixel_svd_scales(
+                self.center, wcs)
+            radius_arcsec = self.radius.to(u.arcsec).value
+            width = 2 * radius_arcsec * scale_major
+            height = 2 * radius_arcsec * scale_minor
+            return EllipsePixelRegion(center, width, height,
+                                      angle=angle,
+                                      meta=self.meta.copy(),
+                                      visual=self.visual.copy())
+
         center, mean_scale = sky_to_pixel_mean_scale(self.center, wcs)
         radius = self.radius.to(u.arcsec).value * mean_scale
         return CirclePixelRegion(center, radius, meta=self.meta.copy(),