Merge pull request #872 from DanRyanIrish/issue871

Fix bug in NDCube.crop

Author: DanRyanIrish

changelog/872.bugfix.rst

@@ -0,0 +1 @@
+Fix bug in `ndcube.NDCube.crop` revealed by trying to crop a 1-D cube. It was caused by the code creating a `~astropy.wcs.wcsapi.SlicedLowLevelWCS` object with a ``slice(None)`` slice item.

ndcube/tests/test_ndcube_slice_and_crop.py

@@ -526,3 +526,21 @@ def test_crop_rotated_celestial(ndcube_4d_ln_lt_l_t):
     small = cube.crop(bottom_left, bottom_right, top_left, top_right)
 
     assert small.data.shape == (1652, 1652)
+
+
+def test_crop_1d():
+    # This use case revealed a bug so has been added as a test.
+    # Create NDCube.
+    wcs = astropy.wcs.WCS(naxis=1)
+    wcs.wcs.ctype = 'WAVE',
+    wcs.wcs.cunit = 'nm',
+    wcs.wcs.cdelt = 4,
+    wcs.wcs.crpix = 1,
+    wcs.wcs.crval = 3,
+    cube = NDCube(np.arange(200), wcs=wcs)
+
+    expected = cube[1:4]
+
+    output = cube.crop((7*u.nm,), (15*u.nm,))
+
+    helpers.assert_cubes_equal(output, expected)

ndcube/utils/cube.py

@@ -141,17 +141,17 @@ def get_crop_item_from_points(points, wcs, crop_by_values, keepdims):
     # where each inner list gives the index of all points for that array axis.
     combined_points_array_idx = [[]] * wcs.pixel_n_dim
     high_level_wcs = HighLevelWCSWrapper(wcs) if isinstance(wcs, BaseLowLevelWCS) else wcs
-    wcs = high_level_wcs.low_level_wcs
+    low_level_wcs = high_level_wcs.low_level_wcs
     # For each point compute the corresponding array indices.
     for point in points:
         # Get the arrays axes associated with each element in point.
         if crop_by_values:
             point_inputs_array_axes = []
-            for i in range(wcs.world_n_dim):
+            for i in range(low_level_wcs.world_n_dim):
                 pix_axes = np.array(
-                    wcs_utils.world_axis_to_pixel_axes(i, wcs.axis_correlation_matrix))
+                    wcs_utils.world_axis_to_pixel_axes(i, low_level_wcs.axis_correlation_matrix))
                 point_inputs_array_axes.append(tuple(
-                    wcs_utils.convert_between_array_and_pixel_axes(pix_axes, wcs.pixel_n_dim)))
+                    wcs_utils.convert_between_array_and_pixel_axes(pix_axes, low_level_wcs.pixel_n_dim)))
             point_inputs_array_axes = tuple(point_inputs_array_axes)
         else:
             point_inputs_array_axes = wcs_utils.array_indices_for_world_objects(high_level_wcs)
@@ -164,14 +164,17 @@ def get_crop_item_from_points(points, wcs, crop_by_values, keepdims):
                 point_indices_with_inputs.append(i)
                 array_axes_with_input.append(point_inputs_array_axes[i])
         array_axes_with_input = set(chain.from_iterable(array_axes_with_input))
-        array_axes_without_input = set(range(wcs.pixel_n_dim)) - array_axes_with_input
+        array_axes_without_input = set(range(low_level_wcs.pixel_n_dim)) - array_axes_with_input
         # Slice out the axes that do not correspond to a coord
         # from the WCS and the input point.
-        wcs_slice = np.array([slice(None)] * wcs.pixel_n_dim)
-        if len(array_axes_without_input):
+        if len(array_axes_without_input) > 0:
+            wcs_slice = np.array([slice(None)] * low_level_wcs.pixel_n_dim)
             wcs_slice[np.array(list(array_axes_without_input))] = 0
-        sliced_wcs = SlicedLowLevelWCS(wcs, slices=tuple(wcs_slice))
-        sliced_point = np.array(point, dtype=object)[np.array(point_indices_with_inputs)]
+            sliced_wcs = SlicedLowLevelWCS(low_level_wcs, slices=tuple(wcs_slice))
+            sliced_point = np.array(point, dtype=object)[np.array(point_indices_with_inputs)]
+        else:
+            # Else, if all axes have at least one crop input, no need to slice the WCS.
+            sliced_wcs, sliced_point = low_level_wcs, np.array(point, dtype=object)
         # Derive the array indices of the input point and place each index
         # in the list corresponding to its axis.
         if crop_by_values: