swev-id: astropy__astropy-13579 – SlicedLowLevelWCS: fix world_to_pixel for sliced axes with coupled PC (fixes #68)

astropy/wcs/setup_package.py

@@ -9,7 +9,13 @@
 from collections import defaultdict
 
 from setuptools import Extension
-from setuptools.dep_util import newer_group
+try:
+    from setuptools.dep_util import newer_group
+except Exception:
+    try:
+        from setuptools._distutils.dep_util import newer_group
+    except Exception:
+        from distutils.dep_util import newer_group
 
 import numpy
 

astropy/wcs/wcsapi/wrappers/sliced_wcs.py

@@ -244,17 +244,39 @@ def pixel_to_world_values(self, *pixel_arrays):
 
     def world_to_pixel_values(self, *world_arrays):
         world_arrays = tuple(map(np.asanyarray, world_arrays))
-        world_arrays_new = []
-        iworld_curr = -1
-        for iworld in range(self._wcs.world_n_dim):
-            if iworld in self._world_keep:
-                iworld_curr += 1
-                world_arrays_new.append(world_arrays[iworld_curr])
-            else:
-                world_arrays_new.append(1.)
 
-        world_arrays_new = np.broadcast_arrays(*world_arrays_new)
-        pixel_arrays = list(self._wcs.world_to_pixel_values(*world_arrays_new))
+        # Fast path: if no world axes were dropped, delegate directly
+        if len(self._world_keep) == self._wcs.world_n_dim:
+            pixel_arrays = list(self._wcs.world_to_pixel_values(*world_arrays))
+        else:
+            # Evaluate the world values on the slice hyperplane for the
+            # dropped world axes. This ensures consistency for correlated WCS.
+            pixel_arrays_eval = []
+            for ipix in range(self._wcs.pixel_n_dim):
+                slc = self._slices_pixel[ipix]
+                if isinstance(slc, numbers.Integral):
+                    pixel_arrays_eval.append(slc)
+                else:
+                    # Any scalar works; dropped world axes are uncorrelated
+                    # with kept pixel axes by construction.
+                    pixel_arrays_eval.append(0)
+
+            pixel_arrays_eval = np.broadcast_arrays(*pixel_arrays_eval)
+            world_eval = self._wcs.pixel_to_world_values(*pixel_arrays_eval)
+
+            # Build full set of world arrays for the underlying call.
+            world_arrays_new = []
+            iworld_curr = -1
+            for iworld in range(self._wcs.world_n_dim):
+                if iworld in self._world_keep:
+                    iworld_curr += 1
+                    world_arrays_new.append(world_arrays[iworld_curr])
+                else:
+                    # Use the world value implied by the fixed slice
+                    world_arrays_new.append(np.asanyarray(world_eval[iworld]))
+
+            world_arrays_new = np.broadcast_arrays(*world_arrays_new)
+            pixel_arrays = list(self._wcs.world_to_pixel_values(*world_arrays_new))
 
         for ipixel in range(self._wcs.pixel_n_dim):
             if isinstance(self._slices_pixel[ipixel], slice) and self._slices_pixel[ipixel].start is not None:

astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py

@@ -899,3 +899,67 @@ def test_pixel_to_world_values_different_int_types():
     for int_coord, np64_coord in zip(int_sliced.pixel_to_world_values(*pixel_arrays),
                                      np64_sliced.pixel_to_world_values(*pixel_arrays)):
         assert all(int_coord == np64_coord)
+
+
+# --- New tests for spectral-spatial mixing via PC ---
+
+def build_coupled_wcs_pc_mix():
+    """
+    Build a simple 3D linear WCS with PC mixing between a spatial
+    world axis (0) and the spectral pixel axis (2).
+
+    Mapping (with CRVAL=0 and CRPIX=0, CDELT=1):
+      w0 = p0 + s * p2
+      w1 = p1
+      w2 = p2
+
+    After slicing out pixel axis 2 (p2 fixed), world axis 2 is dropped,
+    but inverse world->pixel still depends on the correct w2 value.
+    """
+    w = WCS(naxis=3)
+    w.wcs.ctype = ['X', 'Y', 'FREQ']
+    w.wcs.crpix = [0.0, 0.0, 0.0]
+    w.wcs.crval = [0.0, 0.0, 0.0]
+    w.wcs.cdelt = [1.0, 1.0, 1.0]
+    s = 0.3
+    w.wcs.pc = [[1.0, 0.0, s],
+                [0.0, 1.0, 0.0],
+                [0.0, 0.0, 1.0]]
+    return w
+
+
+def test_world_to_pixel_values_coupled_axes_sliced_matches_full():
+    wcs = build_coupled_wcs_pc_mix()
+    x0, y0, z0 = 12.3, -4.7, 3  # z0 integer for exact slicing
+    world_full = wcs.pixel_to_world_values(x0, y0, z0)
+
+    sll = SlicedLowLevelWCS(wcs, [slice(None), slice(None), int(z0)])
+
+    xp, yp = sll.world_to_pixel_values(world_full[0], world_full[1])
+    assert_allclose([xp, yp], [x0, y0], atol=1e-6)
+
+
+def test_pixel_to_world_values_coupled_axes_sliced():
+    wcs = build_coupled_wcs_pc_mix()
+    x0, y0, z0 = 7.5, 2.0, 5
+    world_full = wcs.pixel_to_world_values(x0, y0, z0)
+
+    sll = SlicedLowLevelWCS(wcs, [slice(None), slice(None), int(z0)])
+
+    w0, w1 = sll.pixel_to_world_values(x0, y0)
+    assert_allclose([w0, w1], [world_full[0], world_full[1]], atol=1e-6)
+
+
+def test_world_to_pixel_values_coupled_axes_broadcasting():
+    wcs = build_coupled_wcs_pc_mix()
+    z0 = 4
+    sll = SlicedLowLevelWCS(wcs, [slice(None), slice(None), int(z0)])
+
+    # Build arrays of pixel coords and convert to world, then back.
+    X = np.array([[0.0, 1.5, 3.0], [4.0, 5.5, 6.0]])
+    Y = np.array([[2.0, -1.0, 0.0], [1.0, 2.0, 3.0]])
+    W0, W1 = sll.pixel_to_world_values(X, Y)
+
+    Xr, Yr = sll.world_to_pixel_values(W0, W1)
+    assert_allclose(Xr, X, atol=1e-6)
+    assert_allclose(Yr, Y, atol=1e-6)

docs/changes/wcs/astropy__astropy-13579.bugfix.rst

@@ -0,0 +1,2 @@
+Fix SlicedLowLevelWCS.world_to_pixel_values for sliced correlated WCS (PC mixing): fill dropped world axes with the world values implied by the fixed pixel slice instead of a constant. This corrects large errors in pixel coordinates after slicing.
+