Add some documentation about how crop is inclusive in terms of array indices

[hayesla]

Describe the feature

Following sunpy/sunpy#8276, we should add some defining boundaries in pixel space on which ones are included/not. i.e. the fact that a point on the boundary, it includes that pixel in crop, rather than exclude it, like submap currently does.

Proposed solution

No response

[DanRyanIrish]

Current Rules for NDCube.crop (unless there's a bug)

• Array elements represent an interval in world space in each dimension that is closed on the lower end (includes the lower boundary) and open at the upper end (excludes the boundary).
  • This definition is in keeping with how Python indexing is defined. (Note this is not the same as how WCS defines integer pixel coordinates to be in the centre of the pixel.)
  • This definition ensures that any world coordinate only belongs within one array element, and enables robust, unambiguous cropping behaviour to be implemented.
• NDCube.crop slices the cube to the smallest bounding box in array index space such that all locations in the cube defined by the input world coordinate points are included in the returned cube. This makes NDCube.crop a closed-closed interval operation, i.e. the lower and upper-most array elements are retained.
  • Note that this is same behaviour outlined in the docstring of sunpy.map.GenericMap.submap.
  • Note that this is different to slicing where all array elements from and including the lower array index and up to but not including the upper array index are kept, i.e. it is a closed-open interval operation.
  • The justifications for making crop a closed-closed operation are:
    • it is usually not possible to guarantee that the upper world coordinate will be the open boundary. This is because array elements have a fixed non-zero width in all world dimensions. If we tried to make NDCube.crop closed-closed, and the provided upper world coordinate was somewhere in the middle of an array element, then the actual upper boundary in world space would be an arbitrary distance below the provided world coordinate depending on where it lands in the pixel. By contrast, the closed-open logic of slicing always returns a predictable shape and extent.
    • In the case where both the lower and upper boundaries fall within a single array element, e.g. the 1st, a closed-open logic would mean that the cube would have to both include and exclude that array element. This could result in NDCube trying to crop itself to a scalar, which is not allowed.
    • If the world coordinate was rounded to the nearest array element, we would end up with a behaviour whereby sometimes the provided world coordinate would be included in the cropped cube, and other times not. And it may be hard to predict which way it will go ahead of time with inspected NDCube.wcs.world_to_array_index before the operation.
  • Hence guaranteeing that all locations provided by the user will be in the cropped cube, is a more intuitive and predictable behaviour for users.
  • Note that this behaviour is the

There was confusion regarding this behaviour for NDCube.crop during the 2025 sunpy coordination meeting. However, it was deemed the above rules were true, and the confusion was caused by a bug outlined in #860.

[DanRyanIrish]

The above is describe in the ndcube docs, here. However, it is not explicit and could do with some wordsmithing to be clearer.

[ayshih]

Here's a simple cube to show how NDCube cropping behaves, which I do not feel is the preferred behavior:

import numpy as np

import astropy.units as u
from astropy.wcs import WCS

from ndcube import NDCube

wcs = WCS({
    'ctype1': 'WAVE',
    'cunit1': 'm',
    'cdelt1': 10,
    'crpix1': 1,
    'crval1': 15,
    'ctype2': 'WAVE',
    'cunit2': 'm',
    'cdelt2': 10,
    'crpix2': 1,
    'crval2': 105,
})

cube = NDCube(np.arange(10).reshape(1, 10), wcs=wcs)
cube.crop((40*u.m, 105*u.m), (60*u.m, 105*u.m), keepdims=True)

The cube is defined so that the fourth pixel is 40–50 m and the fifth pixel is 50–60 m. However, since NDCube considers 60 m to "belong" to the sixth pixel (60–70 m), the crop returns three pixels (40–70 m):

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

With PR #874, which fixes the cropping so that pixel edges are not treated as belonging to pixels, the crop returns the desired two pixels (40–60 m):

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

(Incidentally, this cube is specifically constructed so that the floating-point comparisons are exact. It can be difficult in actual cubes for world points to fall exactly on pixel edges.)