astropy · astrofrog · Sep 22, 2026 · Sep 23, 2026 · Sep 23, 2026 · Sep 23, 2026
@@ -17,6 +17,7 @@
 from astropy.utils import unbroadcast
 
 from .wcs import WCS, WCSSUB_LATITUDE, WCSSUB_LONGITUDE
+from .wcsapi.utils import _split_matrix
 
 __doctest_skip__ = ["wcs_to_celestial_frame", "celestial_frame_to_wcs"]
 
@@ -845,45 +846,6 @@ def _pixel_to_pixel_correlation_matrix(wcs_in, wcs_out):
     return matrix
 
 
-def _split_matrix(matrix):
-    """
-    Given an axis correlation matrix from a WCS object, return information about
-    the individual WCS that can be split out.
-
-    The output is a list of tuples, where each tuple contains a list of
-    pixel dimensions and a list of world dimensions that can be extracted to
-    form a new WCS. For example, in the case of a spectral cube with the first
-    two world coordinates being the celestial coordinates and the third
-    coordinate being an uncorrelated spectral axis, the matrix would look like::
-
-        array([[ True,  True, False],
-               [ True,  True, False],
-               [False, False,  True]])
-
-    and this function will return ``[([0, 1], [0, 1]), ([2], [2])]``.
-    """
-    pixel_used = []
-
-    split_info = []
-
-    for ipix in range(matrix.shape[1]):
-        if ipix in pixel_used:
-            continue
-        pixel_include = np.zeros(matrix.shape[1], dtype=bool)
-        pixel_include[ipix] = True
-        n_pix_prev, n_pix = 0, 1
-        while n_pix > n_pix_prev:
-            world_include = matrix[:, pixel_include].any(axis=1)
-            pixel_include = matrix[world_include, :].any(axis=0)
-            n_pix_prev, n_pix = n_pix, np.sum(pixel_include)
-        pixel_indices = list(np.nonzero(pixel_include)[0])
-        world_indices = list(np.nonzero(world_include)[0])
-        pixel_used.extend(pixel_indices)
-        split_info.append((pixel_indices, world_indices))
-
-    return split_info
-
-
 def pixel_to_pixel(wcs_in, wcs_out, *inputs):
     """
     Transform pixel coordinates in a dataset with a WCS to pixel coordinates

@@ -324,12 +324,31 @@ def axis_correlation_matrix(self):
         # each celestial coordinate we copy over the pixel dependencies from
         # the other celestial coordinates.
         celestial = (self.wcs.axis_types // 1000) % 10 == 2
-        celestial_indices = np.nonzero(celestial)[0]
-        for world1 in celestial_indices:
-            for world2 in celestial_indices:
-                if world1 != world2:
-                    matrix[world1] |= matrix[world2]
-                    matrix[world2] |= matrix[world1]
+        matrix[celestial] = matrix[celestial].any(axis=0)
+
+        return matrix
+
+    @property
+    def reverse_axis_correlation_matrix(self):
+        # As for the forward matrix, if there are any distortions present, we
+        # assume that there may be correlations between all axes.
+        if self.has_distortion:
+            return np.ones((self.pixel_n_dim, self.world_n_dim), dtype=bool)
+
+        # Assuming linear world coordinates along each axis, pixel coordinate i
+        # depends on intermediate world coordinate k if the inverse of the PC
+        # matrix is non-zero at (i, k). The numerical inverse of e.g. a rotation
+        # matrix can contain values that are not exactly zero, so we compare
+        # each element to the largest element in its column (all elements in a
+        # column have the same units).
+        inverse_pc = np.abs(np.linalg.inv(self.wcs.get_pc()))
+        matrix = inverse_pc > 1e-10 * inverse_pc.max(axis=0)
+
+        # Each intermediate celestial coordinate is a function of all the
+        # celestial world coordinates, so if a pixel coordinate depends on one
+        # of them it depends on all of them.
+        celestial = (self.wcs.axis_types // 1000) % 10 == 2
+        matrix[:, celestial] = matrix[:, celestial].any(axis=1, keepdims=True)
 
         return matrix
 

@@ -3,6 +3,8 @@
 
 import numpy as np
 
+from .utils import _split_matrix
+
 __all__ = ["BaseLowLevelWCS", "validate_physical_types"]
 
 
@@ -319,6 +321,33 @@ def axis_correlation_matrix(self):
         """
         return np.ones((self.world_n_dim, self.pixel_n_dim), dtype=bool)
 
+    @property
+    def reverse_axis_correlation_matrix(self):
+        """
+        Returns an (`~astropy.wcs.wcsapi.BaseLowLevelWCS.pixel_n_dim`,
+        `~astropy.wcs.wcsapi.BaseLowLevelWCS.world_n_dim`) matrix that
+        indicates using booleans whether a given pixel coordinate depends on a
+        given world coordinate in
+        `~astropy.wcs.wcsapi.BaseLowLevelWCS.world_to_pixel_values`.
+
+        This is not in general the transpose of
+        `~astropy.wcs.wcsapi.BaseLowLevelWCS.axis_correlation_matrix`. By
+        default it is derived from that matrix by treating each pixel
+        coordinate as requiring every world coordinate in the same independent
+        group of axes, which may overstate but never understates the true
+        dependencies. Implementations can override this to return a sparser
+        matrix when fewer world coordinates are needed. Note that when world
+        coordinates carry redundant information there may be several equally
+        valid sparse matrices, and the one returned describes the choice made
+        by the implementation.
+        """
+        # Fill in each independent group of axes, since there is no safe way
+        # to restrict the matrix further without knowing the transformation.
+        reverse = np.zeros((self.pixel_n_dim, self.world_n_dim), dtype=bool)
+        for pixel, world in _split_matrix(self.axis_correlation_matrix):
+            reverse[np.ix_(pixel, world)] = True
+        return reverse
+
     @property
     def serialized_classes(self):
         """

@@ -366,6 +366,18 @@ def test_spectral_cube_nonaligned():
         ],
     )
 
+    # The first pixel axis maps only to latitude, so needs both celestial
+    # world axes but not the frequency, while the other two pixel axes are
+    # mixed with each other by the PC matrix and need all three world axes.
+    assert_equal(
+        wcs.reverse_axis_correlation_matrix,
+        [
+            [True, False, True],
+            [True, True, True],
+            [True, True, True],
+        ],
+    )
+
     # NOTE: we check world_axis_object_components and world_axis_object_classes
     # again here because in the past this failed when non-aligned axes were
     # present, so this serves as a regression test.
@@ -822,28 +834,136 @@ def test_distortion_correlations():
     with pytest.warns(FITSFixedWarning):
         w = WCS(filename)
     assert_equal(w.axis_correlation_matrix, True)
+    assert_equal(w.reverse_axis_correlation_matrix, True)
 
     # Changing PC to an identity matrix doesn't change anything since
     # distortions are still present.
     w.wcs.pc = [[1, 0], [0, 1]]
     assert_equal(w.axis_correlation_matrix, True)
+    assert_equal(w.reverse_axis_correlation_matrix, True)
 
     # Nor does changing the name of the axes to make them non-celestial
     w.wcs.ctype = ["X", "Y"]
     assert_equal(w.axis_correlation_matrix, True)
+    assert_equal(w.reverse_axis_correlation_matrix, True)
 
     # However once we turn off the distortions the matrix changes
     w.sip = None
     assert_equal(w.axis_correlation_matrix, [[True, False], [False, True]])
+    assert_equal(w.reverse_axis_correlation_matrix, [[True, False], [False, True]])
 
     # If we go back to celestial coordinates then the matrix is all True again
     w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
     assert_equal(w.axis_correlation_matrix, True)
+    assert_equal(w.reverse_axis_correlation_matrix, True)
 
     # Or if we change to X/Y but have a non-identity PC
     w.wcs.pc = [[0.9, -0.1], [0.1, 0.9]]
     w.wcs.ctype = ["X", "Y"]
     assert_equal(w.axis_correlation_matrix, True)
+    assert_equal(w.reverse_axis_correlation_matrix, True)
+
+
+# Each tuple gives the CTYPE values, the PC matrix, and the expected reverse
+# axis correlation matrix (as 0/1, converted to bool in the tests).
+# fmt: off
+REVERSE_MATRIX_CASES = [
+    # Independent linear axes: the transpose of the forward matrix
+    (("X", "Y"), [[1, 0], [0, 1]], [[1, 0], [0, 1]]),
+
+    # Lower triangular PC: w0 = p0 and w1 = p0 + p1, so p0 only needs w0
+    (("X", "Y"), [[1, 0], [1, 1]], [[1, 0], [1, 1]]),
+
+    # Upper triangular PC
+    (("X", "Y"), [[1, 1], [0, 1]], [[1, 1], [0, 1]]),
+
+    # Bidiagonal PC whose inverse fills in to a full lower triangle
+    (
+        ("X", "Y", "Z"),
+        [[1, 0, 0], [1, 1, 0], [0, 1, 1]],
+        [[1, 0, 0], [1, 1, 0], [1, 1, 1]],
+    ),
+
+    # Dense lower triangular PC whose inverse is bidiagonal (exact cancellation)
+    (
+        ("X", "Y", "Z"),
+        [[1, 0, 0], [-1, 1, 0], [1, -1, 1]],
+        [[1, 0, 0], [1, 1, 0], [0, 1, 1]],
+    ),
+
+    # Rotated linear axes
+    (("X", "Y"), [[0.9, -0.1], [0.1, 0.9]], [[1, 1], [1, 1]]),
+
+    # Celestial axes always need each other, even with a triangular PC
+    (("RA---TAN", "DEC--TAN"), [[1, 0], [0.3, 1]], [[1, 1], [1, 1]]),
+
+    # Aligned spectral cube
+    (
+        ("RA---TAN", "DEC--TAN", "WAVE"),
+        [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+        [[1, 1, 0], [1, 1, 0], [0, 0, 1]],
+    ),
+
+    # Wavelength skewed by x: the sky pixels never need the wavelength, but
+    # the wavelength pixel needs the sky position to undo the skew
+    (
+        ("RA---TAN", "DEC--TAN", "WAVE"),
+        [[1, 0, 0], [0, 1, 0], [0.1, 0, 1]],
+        [[1, 1, 0], [1, 1, 0], [1, 1, 1]],
+    ),
+
+    # Sky skewed by z: the wavelength pixel never needs the sky position, and
+    # only the skewed sky pixel needs the wavelength to undo the skew
+    (
+        ("RA---TAN", "DEC--TAN", "WAVE"),
+        [[1, 0, 0.1], [0, 1, 0], [0, 0, 1]],
+        [[1, 1, 1], [1, 1, 0], [0, 0, 1]],
+    ),
+
+    # Rastered slit scan: time advances with the first pixel axis, so the time
+    # pixel needs the sky position, but the sky pixels never need the time
+    (
+        ("RA---TAN", "DEC--TAN", "WAVE", "UTC"),
+        [[1, 0, 0, 0], [0.2, 1, 0, 0], [0, 0, 1, 0], [-5.25, 0, 0, 1]],
+        [[1, 1, 0, 0], [1, 1, 0, 0], [0, 0, 1, 0], [1, 1, 0, 1]],
+    ),
+]
+# fmt: on
+
+
+def _reverse_matrix_wcs(ctype, pc):
+    wcs = WCS(naxis=len(ctype))
+    wcs.wcs.ctype = ctype
+    wcs.wcs.crval = [10, 20, 500, 0][: len(ctype)]
+    wcs.wcs.cdelt = [0.01] * len(ctype)
+    wcs.wcs.pc = pc
+    wcs.wcs.set()
+    return wcs
+
+
+@pytest.mark.parametrize(("ctype", "pc", "expected"), REVERSE_MATRIX_CASES)
+def test_reverse_axis_correlation_matrix(ctype, pc, expected):
+    wcs = _reverse_matrix_wcs(ctype, pc)
+    reverse = wcs.reverse_axis_correlation_matrix
+    assert reverse.dtype == bool
+    assert reverse.shape == (wcs.pixel_n_dim, wcs.world_n_dim)
+    assert_equal(reverse, np.array(expected, dtype=bool))
+
+
+@pytest.mark.parametrize(("ctype", "pc"), [case[:2] for case in REVERSE_MATRIX_CASES])
+def test_reverse_axis_correlation_matrix_numerical(ctype, pc):
+    # Perturb each world coordinate in turn (away from the reference point,
+    # where projections are locally diagonal) and check that every pixel
+    # coordinate that responds is marked as depending on that world coordinate.
+    wcs = _reverse_matrix_wcs(ctype, pc)
+    world = np.array(wcs.pixel_to_world_values(*[300, 400, 500, 600][: len(ctype)]))
+    pixel = np.array(wcs.world_to_pixel_values(*world))
+    # Column j of perturbed is the world position with coordinate j perturbed
+    perturbed = world[:, None] + 1e-3 * np.eye(wcs.world_n_dim)
+    changed = ~np.isclose(
+        wcs.world_to_pixel_values(*perturbed), pixel[:, None], rtol=0, atol=1e-9
+    )
+    assert not np.any(changed & ~wcs.reverse_axis_correlation_matrix)
 
 
 def test_custom_ctype_to_ucd_mappings():

@@ -1,6 +1,8 @@
+import numpy as np
 import pytest
+from numpy.testing import assert_equal
 
-from astropy.wcs.wcsapi.low_level_api import validate_physical_types
+from astropy.wcs.wcsapi.low_level_api import BaseLowLevelWCS, validate_physical_types
 
 
 def test_validate_physical_types():
@@ -20,3 +22,83 @@ def test_validate_physical_types():
         ValueError, match=r"'spam' is not a valid IOVA UCD1\+ physical type"
     ):
         validate_physical_types(["spam"])
+
+
+class MatrixLowLevelWCS(BaseLowLevelWCS):
+    """
+    Minimal low-level WCS whose only meaningful property is the axis
+    correlation matrix, used to test the default reverse matrix.
+    """
+
+    def __init__(self, matrix):
+        self._matrix = np.asarray(matrix, dtype=bool)
+
+    @property
+    def pixel_n_dim(self):
+        return self._matrix.shape[1]
+
+    @property
+    def world_n_dim(self):
+        return self._matrix.shape[0]
+
+    @property
+    def axis_correlation_matrix(self):
+        return self._matrix
+
+    @property
+    def world_axis_physical_types(self):
+        return [None] * self.world_n_dim
+
+    @property
+    def world_axis_units(self):
+        return [""] * self.world_n_dim
+
+    # Abstract members that are not needed for these tests
+    pixel_to_world_values = world_to_pixel_values = None
+    world_axis_object_components = world_axis_object_classes = None
+
+
+# Each tuple gives the forward axis correlation matrix and the expected default
+# reverse matrix (as 0/1, converted to bool in the test).
+@pytest.mark.parametrize(
+    ("forward", "expected"),
+    [
+        # Independent axes: the transpose
+        ([[1, 0], [0, 1]], [[1, 0], [0, 1]]),
+        # Fully coupled celestial pair
+        ([[1, 1], [1, 1]], [[1, 1], [1, 1]]),
+        # Spectral cube: two independent blocks
+        ([[1, 1, 0], [1, 1, 0], [0, 0, 1]], [[1, 1, 0], [1, 1, 0], [0, 0, 1]]),
+        # Same cube with the world axes in a different order
+        ([[0, 0, 1], [1, 1, 0], [1, 1, 0]], [[0, 1, 1], [0, 1, 1], [1, 0, 0]]),
+        # Blocks scattered across non-contiguous rows and columns
+        ([[0, 1, 1], [1, 0, 0], [0, 1, 1]], [[0, 1, 0], [1, 0, 1], [1, 0, 1]]),
+        # Triangular: w0 = f(p0), w1 = g(p0, p1), so p1 needs w0 as well as w1
+        ([[1, 0], [1, 1]], [[1, 1], [1, 1]]),
+        # Two pixel, three world axes with a triangular structure
+        ([[1, 0], [1, 1], [1, 1]], [[1, 1, 1], [1, 1, 1]]),
+        # Two pixel, three world axes, all coupled
+        ([[1, 1], [1, 1], [1, 1]], [[1, 1, 1], [1, 1, 1]]),
+        # Triangular block interleaved with an independent axis
+        ([[1, 1, 0], [0, 0, 1], [0, 1, 0]], [[1, 0, 1], [1, 0, 1], [0, 1, 0]]),
+        # World axis that depends on no pixel axis is never required
+        ([[1, 0], [0, 1], [0, 0]], [[1, 0, 0], [0, 1, 0]]),
+        # Pixel axis that no world axis depends on requires nothing
+        ([[1, 0, 0], [0, 1, 0]], [[1, 0], [0, 1], [0, 0]]),
+    ],
+)
+def test_default_reverse_axis_correlation_matrix(forward, expected):
+    wcs = MatrixLowLevelWCS(forward)
+    reverse = wcs.reverse_axis_correlation_matrix
+    assert reverse.dtype == bool
+    assert reverse.shape == (wcs.pixel_n_dim, wcs.world_n_dim)
+    assert_equal(reverse, np.array(expected, dtype=bool))
+
+
+def test_default_reverse_axis_correlation_matrix_all_true():
+    # With no information about the forward matrix, everything is required
+    class AllTrueWCS(MatrixLowLevelWCS):
+        axis_correlation_matrix = BaseLowLevelWCS.axis_correlation_matrix
+
+    wcs = AllTrueWCS(np.ones((3, 2)))
+    assert_equal(wcs.reverse_axis_correlation_matrix, np.ones((2, 3), dtype=bool))