Error when cropping a 1D NDCube

[ayshih]

Describe the bug

Attempting to crop a 1D NDcube errors.

To Reproduce

import numpy as np

import astropy.wcs
import astropy.units as u

from ndcube import NDCube

data_1d = np.arange(200)

wcs_1d = astropy.wcs.WCS(naxis=1)
wcs_1d.wcs.ctype = 'WAVE',
wcs_1d.wcs.cunit = 'nm',
wcs_1d.wcs.cdelt = 4,
wcs_1d.wcs.crpix = 1,
wcs_1d.wcs.crval = 3,

cube_1d = NDCube(data_1d, wcs=wcs_1d)
cube_1d.crop((7*u.nm,), (15*u.nm,))

gives

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[78], line 1
----> 1 cube_1d.crop((7*u.nm,), (15*u.nm,))

File ~\AppData\Local\mambaforge\envs\py312\Lib\site-packages\ndcube\ndcube.py:622, in NDCubeBase.crop(self, wcs, keepdims, *points)
    619 def crop(self, *points, wcs=None, keepdims=False):
    620     # The docstring is defined in NDCubeABC
    621     # Calculate the array slice item corresponding to bounding box and return sliced cube.
--> 622     item = self._get_crop_item(*points, wcs=wcs, keepdims=keepdims)
    623     return self[item]

File ~\AppData\Local\mambaforge\envs\py312\Lib\site-packages\ndcube\utils\cube.py:57, in sanitize_wcs.<locals>.wcs_wrapper(*args, **kwargs)
     53     raise TypeError("wcs argument must be a High Level WCS or an ExtraCoords object.")
     55 params.arguments['wcs'] = wcs
---> 57 return func(*params.args, **params.kwargs)

File ~\AppData\Local\mambaforge\envs\py312\Lib\site-packages\ndcube\ndcube.py:648, in NDCubeBase._get_crop_item(self, wcs, keepdims, *points)
    644         if not (value is None or isinstance(value, classes[j])):
    645             raise TypeError(f"{type(value)} of component {j} in point {i} is "
    646                             f"incompatible with WCS component {comp[j]} "
    647                             f"{classes[j]}.")
--> 648 return utils.cube.get_crop_item_from_points(points, wcs, False, keepdims=keepdims)

File ~\AppData\Local\mambaforge\envs\py312\Lib\site-packages\ndcube\utils\cube.py:183, in get_crop_item_from_points(points, wcs, crop_by_values, keepdims)
    181         point_array_indices = tuple(a.item() for a in point_array_indices)
    182 else:
--> 183     point_array_indices = HighLevelWCSWrapper(sliced_wcs).world_to_array_index(
    184         *sliced_point)
    185     # If returned value is a 0-d array, convert to a length-1 tuple.
    186     if isinstance(point_array_indices, np.ndarray) and point_array_indices.ndim == 0:

File ~\OneDrive - NASA\code\astropy\astropy\wcs\wcsapi\high_level_api.py:129, in BaseHighLevelWCS.world_to_array_index(self, *world_objects)
    117 """
    118 Convert world coordinates (represented by Astropy objects) to array
    119 indices.
   (...)
    126 as rounded integers.
    127 """
    128 if self.low_level_wcs.pixel_n_dim == 1:
--> 129     return _toindex(self.world_to_pixel(*world_objects))
    130 else:
    131     return tuple(_toindex(self.world_to_pixel(*world_objects)[::-1]).tolist())

File ~\OneDrive - NASA\code\astropy\astropy\wcs\wcsapi\high_level_api.py:356, in HighLevelWCSMixin.world_to_pixel(self, *world_objects)
    351 world_values = high_level_objects_to_values(
    352     *world_objects, low_level_wcs=self.low_level_wcs
    353 )
    355 # Finally we convert to pixel coordinates
--> 356 pixel_values = self.low_level_wcs.world_to_pixel_values(*world_values)
    358 return pixel_values

File ~\OneDrive - NASA\code\astropy\astropy\wcs\wcsapi\wrappers\sliced_wcs.py:271, in SlicedLowLevelWCS.world_to_pixel_values(self, *world_arrays)
    268         world_arrays_new.append(sliced_out_world_coords[iworld])
    270 world_arrays_new = np.broadcast_arrays(*world_arrays_new)
--> 271 pixel_arrays = list(self._wcs.world_to_pixel_values(*world_arrays_new))
    273 for ipixel in range(self._wcs.pixel_n_dim):
    274     if (
    275         isinstance(self._slices_pixel[ipixel], slice)
    276         and self._slices_pixel[ipixel].start is not None
    277     ):

TypeError: iteration over a 0-d array

Adding a second dimension to the NDCube "fixes" the problem:

data_2d = data_1d[np.newaxis, :]

wcs_2d = astropy.wcs.WCS(naxis=2)
wcs_2d.wcs.ctype = 'WAVE', 'WAVE'
wcs_2d.wcs.cunit = 'nm', 'nm'
wcs_2d.wcs.cdelt = 4, 1
wcs_2d.wcs.crpix = 1, 1
wcs_2d.wcs.crval = 3, 1

cube_2d = NDCube(data_2d, wcs=wcs_2d)
cube_2d.crop((7*u.nm, 1*u.nm), (15*u.nm, 1*u.nm))

gives

<ndcube.ndcube.NDCube object at 0x00000165E2AE9490>
NDCube
------
Shape: (3,)
Physical Types of Axes: [('em.wl',)]
Unit: None
Data Type: int64

Screenshots

No response

System Details

==============================
sunpy Installation Information

General
#######
OS: Windows 11 10.0.26100
Arch: 64bit, (Intel64 Family 6 Model 85 Stepping 7, GenuineIntel)
sunpy: 6.1.dev119+g7259b5058.d20240911
Installation path: C:\Users\ayshih\AppData\Local\mambaforge\envs\py312\Lib\site-packages\sunpy-6.1.dev119+g7259b5058.d20240911.dist-info

Required Dependencies
#####################
astropy: 7.0.0.dev810+gb982985335.d20240911
numpy: 2.1.1
packaging: 24.1
parfive: 2.1.0
pyerfa: 2.0.1.4
requests: 2.32.3

Optional Dependencies
#####################
sunpy: 6.1.dev119+g7259b5058.d20240911

Installation method

No response

[DanRyanIrish]

This appears to be caused because the code here creates a SlicedLowLevelWCS with slice item that in this case is slice(None, None, None), i.e. no slicing is actually done. That appears to cause a problem when we convert the WCS back to a HighLevelWCSWrapper and call world_to_array_index here.

I'm not sure whether to classify this as an astropy bug. I would have thought this should work fine. Nonetheless, it's something we can fix in the ndcube code base.

[ayshih]

I'm not sure whether to classify this as an astropy bug. I would have thought this should work fine.

It's indeed a bug, and I just put up the fix (astropy/astropy#18394)