Fix crop for WCSes with interleaved world objects

changelog/977.bugfix.1.rst

@@ -0,0 +1 @@
+ ~ndcube.NDCube.crop and ~ndcube.NDCube.crop_by_values no longer swap the bounds of some axes on cubes with nine or more dimensions.

changelog/977.bugfix.rst

@@ -0,0 +1 @@
+~ndcube.NDCube.crop now accepts points in the order ``pixel_to_world`` returns, including partial points with `None`, when a WCS's world objects are interleaved (e.g. HPLN, WAVE, HPLT). If each object matches only one world object, the objects can be given in any order, as with ~astropy.wcs.wcsapi.BaseHighLevelWCS.world_to_pixel. Errors for a point of the wrong length or type now list the expected objects.

ndcube/ndcube.py

@@ -633,22 +633,26 @@ def _get_crop_item(self, *points, wcs=None, keepdims=False):
         # Quit out early if we are no-op
         if no_op:
             return tuple([slice(None)] * wcs.pixel_n_dim)
-        comp = [c[0] for c in wcs.world_axis_object_components]
-        # Trim to unique component names - `np.unique(..., return_index=True)
-        # keeps sorting alphabetically, set() seems just nondeterministic.
-        for k, c in enumerate(comp):
-            if comp.count(c) > 1:
-                comp.pop(k)
+        comp = utils.misc.unique_sorted(c[0] for c in wcs.world_axis_object_components)
         classes = [wcs.world_axis_object_classes[c][0] for c in comp]
+        expected = ", ".join(f"{name} ({cls.__name__})" for name, cls in zip(comp, classes))
         for i, point in enumerate(points):
             if len(point) != len(comp):
                 raise ValueError(f"{len(point)} components in point {i} do not match "
-                                 f"WCS with {len(comp)} components.")
+                                 f"WCS with {len(comp)} components. Each point must "
+                                 "have one entry per world object (use None for a "
+                                 f"component that should not be cropped), in order: {expected}.")
+            # Like astropy's world_to_pixel, match objects to components by class when unambiguous (#608).
+            vals = [v for v in point if v is not None]
+            slots = [[j for j, cls in enumerate(classes) if isinstance(v, cls)] for v in vals]
+            matched = {s[0]: v for v, s in zip(vals, slots) if len(s) == 1}
+            if len(matched) == len(vals):
+                points[i] = point = [matched.get(j) for j in range(len(comp))]
             for j, value in enumerate(point):
                 if not (value is None or isinstance(value, classes[j])):
                     raise TypeError(f"{type(value)} of component {j} in point {i} is "
                                     f"incompatible with WCS component {comp[j]} "
-                                    f"{classes[j]}.")
+                                    f"{classes[j]}. Expected order: {expected}.")
         return utils.cube.get_crop_item_from_points(points, wcs, False, keepdims=keepdims,
                                                     original_shape=self.data.shape)
 

ndcube/tests/test_ndcube_slice_and_crop.py

@@ -287,7 +287,8 @@ def test_crop_missing_dimensions(ndcube_4d_ln_lt_l_t):
     interval0 = cube.wcs.array_index_to_world([1, 2], [0, 1], [0, 1], [0, 2])[0]
     lower_corner = [interval0[0], None]
     upper_corner = [interval0[-1], None]
-    with pytest.raises(ValueError, match=r'2 components in point 0 do not match WCS with 3'):
+    with pytest.raises(ValueError, match=r'2 components in point 0 do not match WCS with 3 .* in order: '
+                                         r'time \(Time\), spectral \(Quantity\), celestial \(SkyCoord\)\.$'):
         cube.crop(lower_corner, upper_corner)
 
 
@@ -299,7 +300,8 @@ def test_crop_mismatch_class(ndcube_4d_ln_lt_l_t):
     lower_corner = [coord[0] for coord in intervals]
     upper_corner = [coord[-1] for coord in intervals]
     with pytest.raises(TypeError, match=r"<class .*.SpectralCoord'> of component 0 in point 0 is "
-                                        r"incompatible with WCS component time"):
+                                        r"incompatible with WCS component time .* Expected order: "
+                                        r"time \(Time\), spectral \(Quantity\), celestial \(SkyCoord\)\.$"):
         cube.crop(lower_corner, upper_corner)
 
 
@@ -470,6 +472,9 @@ def test_crop_by_extra_coords_all_axes_with_coord(ndcube_3d_ln_lt_l_ec_all_axes)
     output = cube.crop(lower_corner, upper_corner, wcs=cube.extra_coords)
     expected = cube[0, 0:2, 1:4]
     helpers.assert_cubes_equal(output, expected)
+    # 1 m matches both Quantity coords, so the point is matched by position.
+    output = cube.crop((None, None, interval2[0]), (None, None, interval2[1]), wcs=cube.extra_coords)
+    helpers.assert_cubes_equal(output, cube[:, :, 1:4])
 
 
 def test_crop_by_extra_coords_values_all_axes_with_coord(ndcube_3d_ln_lt_l_ec_all_axes):
@@ -617,6 +622,36 @@ def test_crop_all_points_beyond_cube_extent_error(points):
         cube.crop(*points, keepdims=True)
 
 
+def test_crop_non_contiguous_world_objects(ndcube_4d_ln_l_t_lt):
+    # World axes are HPLT, TIME, WAVE, HPLN, so the SkyCoord components are not
+    # adjacent. The expected object order is first appearance: SkyCoord, Time, Quantity.
+    cube = ndcube_4d_ln_l_t_lt
+    lower = cube.wcs.array_index_to_world(1, 0, 0, 4)[0]
+    upper = cube.wcs.array_index_to_world(2, 0, 0, 5)[0]
+    output = cube.crop([lower, None, None], [upper, None, None])
+    helpers.assert_cubes_equal(output, cube[1:3, :, :, 4:6])
+    # Unambiguous objects are matched to components by class, so any order works (#608).
+    sky, _, wave = cube.wcs.pixel_to_world(1, 2, 3, 4)
+    helpers.assert_cubes_equal(cube.crop([wave, sky, None], keepdims=True), cube[4:5, 3:4, :, 1:2])
+    # 1 matches no class, so the point is matched by position and rejected.
+    with pytest.raises(TypeError):
+        cube.crop([wave, sky, 1])
+
+
+def test_crop_by_values_many_axes_keeps_axis_order():
+    # Pixel axes 3 and 8 used to be visited in set order ([8, 3]) and swapped.
+    wcs = WCS(naxis=9)
+    wcs.wcs.ctype = ["LINEAR"] * 9
+    wcs.wcs.cunit = ["m"] * 9
+    cube = NDCube(np.zeros((5,) * 9, dtype=bool), wcs=wcs)
+    point = [None] * 9
+    point[3] = 1 * u.m
+    point[8] = 4 * u.m
+    output = cube.crop_by_values(point, keepdims=True)
+    world = output.wcs.low_level_wcs.pixel_to_world_values(*[0] * 9)
+    assert (world[3], world[8]) == (1, 4)
+
+
 def test_crop_by_values_quantity_table_coordinate():
     # Regression: QuantityTableCoordinate-based WCS raised
     # "High Level objects are not supported with the native API" because

ndcube/utils/cube.py

@@ -164,7 +164,7 @@ def get_crop_item_from_points(points, wcs, crop_by_values, keepdims, original_sh
                 pixel_axes_with_input.append(point_inputs_pixel_axes[i])
         pixel_axes_with_input = set(chain.from_iterable(pixel_axes_with_input))
         pixel_axes_without_input = set(range(low_level_wcs.pixel_n_dim)) - pixel_axes_with_input
-        pixel_axes_with_input = np.array(list(pixel_axes_with_input))
+        pixel_axes_with_input = np.array(sorted(pixel_axes_with_input))
         pixel_axes_without_input = np.array(list(pixel_axes_without_input))
         # Slice out the axes that do not correspond to a coord
         # from the WCS and the input point.