[alasdairwilson]

currently machine precision means that we are only correctly finding pixel edges if A) we do an exact binary sum or B) 50/50 chance the machine accuracy lies on the correct side of the boundary.

By complete luck our existing tests for this were either all binary exact (i.e. exact same units in world coord as WCS AND cdelt a power of 2 divisor) or they were lucky to round to the correct side of the pixel.

Failures though are easy to find:

import numpy as np
from astropy.wcs import WCS
from ndcube import NDCube

cube = NDCube(np.arange(120).reshape(10, 12),
              wcs=WCS({"CTYPE1": "RA---TAN", "CTYPE2": "DEC--TAN"}))

# Corners on the pixel edges at 0.5 and 2.5, so this is a 2x2 pixel region.
corners = cube.wcs.pixel_to_world([0.5, 2.5, 2.5, 0.5], [0.5, 0.5, 2.5, 2.5])
x, y = cube.wcs.world_to_pixel(corners)

print("crop    :", cube.crop(*corners, keepdims=True).shape)

returns a 2,4 area instead of 2,2

even more disturbing:

wcs = WCS({"CTYPE1": "WAVE", "CUNIT1": "Angstrom",
           "CDELT1": 0.2, "CRPIX1": 0, "CRVAL1": 0})
cube = NDCube(np.arange(10), wcs=wcs)

gives inconsistent results given the units of the spectralcoord that accesses that axis

cube.crop([0.3 * u.Angstrom], [0.45 * u.Angstrom], keepdims=True)
cube.crop([0.03 * u.nm],      [0.45 * u.Angstrom], keepdims=True)

the first has shape 1 and the second 1,0 as the first rounds to 0.5000000002 and the second to 0.499999998

AI tools were used for:

• No AI tools were used

Regardless of AI use, the human contributor remains fully responsible for correctness, design choices, licensing compatibility, and long-term maintainability.

[ayshih]

I have not reviewed the greater context for this change, but I would discourage the use of np.round() here because it rounds half to even, and I suspect that rounding half up would be better here. (An example of rounding half up is on the very next line: np.floor(... + 0.5).)

[alasdairwilson]

In this case what it is doing is not rounding to even integers, it would round 1.5e-6 to 2e-6 and 2.5e-6 to 2e-6 (rather than 3e-6 that half-up would do).

I'm ok with this behaviour, especially since what are the chances when we are dealing with this machine precision issue at 1e-11 that after the wcs transform we get an exact round value at the 7th decimal: i.e. a number that looks like: 0.xxxxxx50000000. In fact I am pretty sure getting such a value is structurally impossible given we are dealing with machine precision 5+ decimal places smaller.

And a line of code doing floor/ceiling on the 7th decimal of a figure is gonna look insane to anyone following along after (as in doing np.floor(10**7*x+0.5)/10**7 )

[ayshih]

I'm ok with this behaviour, especially since what are the chances when we are dealing with this machine precision issue at 1e-11 that after the wcs transform we get an exact round value at the 7th decimal: i.e. a number that looks like: 0.xxxxxx50000000. In fact I am pretty sure getting such a value is structurally impossible given we are dealing with machine precision 5+ decimal places smaller.

It doesn't have to exact 5 at the 7th decimal, and as a formal matter, such a number can't be stored in a Python float. It just has to be exact 0.5 when multiplied 1e6. It seems straightforward to construct an example of a world coordinate that is effectively exactly 5 at the 7th decimal in pixel coordinates.

My recommendation is to not do this pre-rounding that you are suggesting, but to instead simply modify the following two lines to:

            min_array_idx = int(np.floor(min(pixel_coords) + 0.5 + 5e-7))
            max_array_idx = int(np.ceil(max(pixel_coords) - 0.5 - 5e-7)) + 1

This is the approach I have used in matplotlib for similar situations (except I happen to use 1e-8).

[alasdairwilson]

This is now done