Use principal WCS axes

Author: larrybradley

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

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,
@@ -179,12 +181,12 @@ def to_sky(self, wcs, *, use_ellipse=False):
             ``use_ellipse`` is `True`.
         """
         if use_ellipse:
-            center, scale_w, scale_h, angle = pixel_to_sky_scales(
-                self.center, wcs, 0.0)
-            inner_width = 2 * self.inner_radius * scale_w * u.arcsec
-            outer_width = 2 * self.outer_radius * scale_w * u.arcsec
-            inner_height = 2 * self.inner_radius * scale_h * u.arcsec
-            outer_height = 2 * self.outer_radius * scale_h * u.arcsec
+            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,
@@ -286,14 +288,14 @@ def to_pixel(self, wcs, *, use_ellipse=False):
             ``use_ellipse`` is `True`.
         """
         if use_ellipse:
-            center, scale_w, scale_h, angle = sky_to_pixel_scales(
-                self.center, wcs, 0.0)
+            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_w
-            outer_width = 2 * outer_radius_arcsec * scale_w
-            inner_height = 2 * inner_radius_arcsec * scale_h
-            outer_height = 2 * outer_radius_arcsec * scale_h
+            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,

regions/shapes/circle.py

@@ -11,9 +11,9 @@
 
 from regions._geometry import circular_overlap_grid
 from regions._utils.wcs_helpers import (pixel_to_sky_mean_scale,
-                                        pixel_to_sky_scales,
+                                        pixel_to_sky_svd_scales,
                                         sky_to_pixel_mean_scale,
-                                        sky_to_pixel_scales)
+                                        sky_to_pixel_svd_scales)
 from regions.core.attributes import (PositiveScalar, PositiveScalarAngle,
                                      RegionMetaDescr, RegionVisualDescr,
                                      ScalarPixCoord, ScalarSkyCoord)
@@ -115,10 +115,10 @@ def to_sky(self, wcs, *, use_ellipse=False):
             is `True`.
         """
         if use_ellipse:
-            center, scale_w, scale_h, angle = pixel_to_sky_scales(
-                self.center, wcs, 0.0)
-            width = Angle(2 * self.radius * scale_w, 'arcsec')
-            height = Angle(2 * self.radius * scale_h, 'arcsec')
+            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())
@@ -295,11 +295,11 @@ def to_pixel(self, wcs, *, use_ellipse=False):
             is `True`.
         """
         if use_ellipse:
-            center, scale_w, scale_h, angle = sky_to_pixel_scales(
-                self.center, wcs, 0.0)
+            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_w
-            height = 2 * radius_arcsec * scale_h
+            width = 2 * radius_arcsec * scale_major
+            height = 2 * radius_arcsec * scale_minor
             return EllipsePixelRegion(center, width, height,
                                       angle=angle,
                                       meta=self.meta.copy(),

regions/shapes/tests/test_circle.py

@@ -169,8 +169,8 @@ def test_to_sky_ellipse_roundtrip(self):
         pix_ellipse = sky_ellipse.to_pixel(self.wcs)
         # For a simple WCS without distortion, the roundtrip should
         # recover the original diameter as width and height
-        assert_allclose(pix_ellipse.width, 2 * self.radius, rtol=1e-5)
-        assert_allclose(pix_ellipse.height, 2 * self.radius, rtol=1e-5)
+        assert_allclose(pix_ellipse.width, 2 * self.radius, rtol=1e-4)
+        assert_allclose(pix_ellipse.height, 2 * self.radius, rtol=1e-4)
         assert_allclose(pix_ellipse.center.x, self.center.x, rtol=1e-5)
         assert_allclose(pix_ellipse.center.y, self.center.y, rtol=1e-5)