diff --git a/astropy/wcs/setup_package.py b/astropy/wcs/setup_package.py index 3b7b416d5eff..1165aff4bce7 100644 --- a/astropy/wcs/setup_package.py +++ b/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 diff --git a/astropy/wcs/wcsapi/wrappers/sliced_wcs.py b/astropy/wcs/wcsapi/wrappers/sliced_wcs.py index d7605b078c31..3bfb2c654906 100644 --- a/astropy/wcs/wcsapi/wrappers/sliced_wcs.py +++ b/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: diff --git a/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py b/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py index a3541149624e..075ddb8adf8e 100644 --- a/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py +++ b/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) diff --git a/docs/changes/wcs/astropy__astropy-13579.bugfix.rst b/docs/changes/wcs/astropy__astropy-13579.bugfix.rst new file mode 100644 index 000000000000..31d84bd2dc67 --- /dev/null +++ b/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. +