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

Activity

1. hayesla

   added

   Feature RequestNew feature wanted.New feature wanted.

   on Jun 19, 2025
2. hayesla

   changed the title [-]Add some documentation about how crop is inclusive in terms of pixels[/-] [+]Add some documentation about how crop is inclusive in terms of pixel coords[/+] on Jun 19, 2025
3. DanRyanIrish

   changed the title [-]Add some documentation about how crop is inclusive in terms of pixel coords[/-] [+]Add some documentation about how crop is inclusive in terms of array indices[/+] on Jun 30, 2025
4. DanRyanIrish commented on Jun 30, 2025

   

   DanRyanIrish

   on Jun 30, 2025

   Last edited by DanRyanIrish

   Member

   More actions

   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.
5. DanRyanIrish commented on Jun 30, 2025

   

   DanRyanIrish

   on Jun 30, 2025

   Member

   More actions

   The above is describe in the ndcube docs, here. However, it is not explicit and could do with some wordsmithing to be clearer.
6. ayshih commented on Jul 8, 2025

   

   ayshih

   on Jul 8, 2025

   Last edited by ayshih

   Member

   More actions

   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.)