Repository: astropy/astropy

Issue:
Cutout2D should include a parameter for the limit rounding method
### What is the problem this feature will solve?

When making cutouts with `Cutout2D`, the pixel limits are calculated using `np.ceil`. For a small cutout, this could potentially create a large amount of offset from the requested center. For example, if the unrounded pixel limit were calculated to be 2.00001, the better choice would be to round down to 2 rather than up to 3. It would be great if `Cutout2D` had a parameter, something like `limit_rounding_method`, that could default to `np.ceil` but also be set to `round` for `np.round` and `floor` for `np.floor`. 

### Describe the desired outcome

The cutout's pixel limits are rounded according to the value of the parameter. In the case that a cutout is so small that it rounds to 0, the pixel limits should be set to include the pixel that the requested coordinate falls within.

### Additional context

_No response_

Discussion:

Patch:
diff --git a/astropy/nddata/utils.py b/astropy/nddata/utils.py
index 0474bcc8fc..a2c9ef4cb5 100644
--- a/astropy/nddata/utils.py
+++ b/astropy/nddata/utils.py
@@ -34,7 +34,14 @@ class PartialOverlapError(ValueError):
     """Raised when arrays only partially overlap."""
 
 
-def overlap_slices(large_array_shape, small_array_shape, position, mode="partial"):
+def overlap_slices(
+    large_array_shape,
+    small_array_shape,
+    position,
+    mode="partial",
+    *,
+    limit_rounding_method=np.ceil,
+):
     """
     Get slices for the overlapping part of a small and a large array.
 
@@ -69,6 +76,10 @@ def overlap_slices(large_array_shape, small_array_shape, position, mode="partial
         otherwise an `~astropy.nddata.utils.PartialOverlapError` is
         raised.  In all modes, non-overlapping arrays will raise a
         `~astropy.nddata.utils.NoOverlapError`.
+    limit_rounding_method : callable
+        The rounding method when calculating the minimum and maximum pixel indices.
+        This must be a callable function. Examples: `~numpy.ceil`, `~numpy.floor`,
+        `~numpy.round`. Default is `~numpy.ceil`.
 
     Returns
     -------
@@ -104,14 +115,23 @@ def overlap_slices(large_array_shape, small_array_shape, position, mode="partial
             '"position" must have the same number of dimensions as "small_array_shape".'
         )
 
+    if not callable(limit_rounding_method):
+        raise ValueError("Limit rounding method must be a callable function.")
+
     # define the min/max pixel indices
-    indices_min = [
-        int(np.ceil(pos - (small_shape / 2.0)))
-        for (pos, small_shape) in zip(position, small_array_shape)
-    ]
+    # round according to the limit_rounding_method
+    try:
+        indices_min = [
+            int(limit_rounding_method(pos - (small_shape / 2.0)))
+            for (pos, small_shape) in zip(position, small_array_shape)
+        ]
+    except (TypeError, ValueError) as exc:
+        raise ValueError(
+            "Limit rounding method must accept a single number as input and return a single number."
+        ) from exc
     indices_max = [
-        int(np.ceil(pos + (small_shape / 2.0)))
-        for (pos, small_shape) in zip(position, small_array_shape)
+        int(idx_min + small_shape)
+        for (idx_min, small_shape) in zip(indices_min, small_array_shape)
     ]
 
     for e_max in indices_max:
@@ -159,6 +179,8 @@ def extract_array(
     mode="partial",
     fill_value=np.nan,
     return_position=False,
+    *,
+    limit_rounding_method=np.ceil,
 ):
     """
     Extract a smaller array of the given shape and position from a
@@ -199,6 +221,10 @@ def extract_array(
     return_position : bool, optional
         If `True`, return the coordinates of ``position`` in the
         coordinate system of the returned array.
+    limit_rounding_method : callable
+        The rounding method when calculating the minimum and maximum pixel indices.
+        This must be a callable function. Examples: `~numpy.ceil`, `~numpy.floor`,
+        `~numpy.round`. Default is `~numpy.ceil`.
 
     Returns
     -------
@@ -234,7 +260,11 @@ def extract_array(
         raise ValueError("Valid modes are 'partial', 'trim', and 'strict'.")
 
     large_slices, small_slices = overlap_slices(
-        array_large.shape, shape, position, mode=mode
+        array_large.shape,
+        shape,
+        position,
+        mode=mode,
+        limit_rounding_method=limit_rounding_method,
     )
     extracted_array = array_large[large_slices]
     if return_position:
@@ -264,7 +294,7 @@ def extract_array(
         return extracted_array
 
 
-def add_array(array_large, array_small, position):
+def add_array(array_large, array_small, position, *, limit_rounding_method=np.ceil):
     """
     Add a smaller array at a given position in a larger array.
 
@@ -278,6 +308,10 @@ def add_array(array_large, array_small, position):
     position : tuple
         Position of the small array's center, with respect to the large array.
         Coordinates should be in the same order as the array shape.
+    limit_rounding_method : callable
+        The rounding method when calculating the minimum and maximum pixel indices.
+        This must be a callable function. Examples: `~numpy.ceil`, `~numpy.floor`,
+        `~numpy.round`. Default is `~numpy.ceil`.
 
     Returns
     -------
@@ -311,7 +345,10 @@ def add_array(array_large, array_small, position):
         for (large_shape, small_shape) in zip(array_large.shape, array_small.shape)
     ):
         large_slices, small_slices = overlap_slices(
-            array_large.shape, array_small.shape, position
+            array_large.shape,
+            array_small.shape,
+            position,
+            limit_rounding_method=limit_rounding_method,
         )
         array_large[large_slices] += array_small[small_slices]
         return array_large
@@ -448,6 +485,11 @@ class Cutout2D:
         into the original ``data`` array.  If `True`, then the
         cutout data will hold a copy of the original ``data`` array.
 
+    limit_rounding_method : callable
+        The rounding method when calculating the minimum and maximum pixel indices.
+        This must be a callable function. Examples: `~numpy.ceil`, `~numpy.floor`,
+        `~numpy.round`. Default is `~numpy.ceil`.
+
     Attributes
     ----------
     data : 2D `~numpy.ndarray`
@@ -541,7 +583,16 @@ class Cutout2D:
     """
 
     def __init__(
-        self, data, position, size, wcs=None, mode="trim", fill_value=np.nan, copy=False
+        self,
+        data,
+        position,
+        size,
+        wcs=None,
+        mode="trim",
+        fill_value=np.nan,
+        copy=False,
+        *,
+        limit_rounding_method=np.ceil,
     ):
         if wcs is None:
             wcs = getattr(data, "wcs", None)
@@ -603,6 +654,7 @@ class Cutout2D:
             mode=mode,
             fill_value=fill_value,
             return_position=True,
+            limit_rounding_method=limit_rounding_method,
         )
         if copy:
             cutout_data = np.copy(cutout_data)
@@ -610,7 +662,11 @@ class Cutout2D:
 
         self.input_position_cutout = input_position_cutout[::-1]  # (x, y)
         slices_original, slices_cutout = overlap_slices(
-            data.shape, shape, pos_yx, mode=mode
+            data.shape,
+            shape,
+            pos_yx,
+            mode=mode,
+            limit_rounding_method=limit_rounding_method,
         )
 
         self.slices_original = slices_original
diff --git a/docs/changes/nddata/17876.feature.rst b/docs/changes/nddata/17876.feature.rst
new file mode 100644
index 0000000000..f48f9f9844
--- /dev/null
+++ b/docs/changes/nddata/17876.feature.rst
@@ -0,0 +1,5 @@
+Add the ``limit_rounding_method`` parameter to `~astropy.nddata.Cutout2D`, 
+`~astropy.nddata.overlap_slices`, `~astropy.nddata.extract_array`, and 
+`~astropy.nddata.add_array` to allow users to specify the rounding method 
+used when calculating the pixel limits of the cutout. The default method 
+is to use `~numpy.ceil`.
\ No newline at end of file


Test patch:
diff --git a/astropy/nddata/tests/test_utils.py b/astropy/nddata/tests/test_utils.py
index c160177ba3..ce41283f2f 100644
--- a/astropy/nddata/tests/test_utils.py
+++ b/astropy/nddata/tests/test_utils.py
@@ -148,6 +148,84 @@ def test_slices_nonfinite_position(position):
         overlap_slices((7, 7), (3, 3), position)
 
 
+@pytest.mark.parametrize(
+    "small_array_shape, position, limit_rounding_method, expected_lg, expected_sm",
+    [
+        (
+            (3, 5),
+            (4, 5),
+            np.ceil,
+            (slice(3, 6), slice(3, 8)),
+            (slice(0, 3), slice(0, 5)),
+        ),
+        (
+            (3, 5),
+            (4, 5),
+            np.floor,
+            (slice(2, 5), slice(2, 7)),
+            (slice(0, 3), slice(0, 5)),
+        ),
+        (
+            (3, 5),
+            (4.2, 5.6),
+            np.round,
+            (slice(3, 6), slice(3, 8)),
+            (slice(0, 3), slice(0, 5)),
+        ),
+        (
+            (3, 3),
+            (6, 6),
+            np.round,
+            (slice(4, 7), slice(4, 7)),
+            (slice(0, 3), slice(0, 3)),
+        ),
+        (
+            (3, 3),
+            (5, 5),
+            np.round,
+            (slice(4, 7), slice(4, 7)),
+            (slice(0, 3), slice(0, 3)),
+        ),
+        (
+            (3, 3),
+            (5, 5),
+            np.trunc,
+            (slice(3, 6), slice(3, 6)),
+            (slice(0, 3), slice(0, 3)),
+        ),
+    ],
+)
+def test_slices_limit_rounding_method(
+    small_array_shape, position, limit_rounding_method, expected_lg, expected_sm
+):
+    """Call overlap_slices with different limit rounding methods."""
+    slc_lg, slc_sm = overlap_slices(
+        (10, 10),
+        small_array_shape,
+        position,
+        limit_rounding_method=limit_rounding_method,
+    )
+    assert slc_lg == expected_lg
+    assert slc_sm == expected_sm
+
+
+@pytest.mark.parametrize(
+    "limit_rounding_method, error_msg",
+    [
+        (None, "^Limit rounding method must be a callable"),
+        ("invalid", "^Limit rounding method must be a callable"),
+        (lambda x: (1, 2), "^Limit rounding method must accept a single number"),
+        (lambda x: "invalid", "^Limit rounding method must accept a single number"),
+        (lambda x, y: 1, "^Limit rounding method must accept a single number"),
+        (lambda: 1, "^Limit rounding method must accept a single number"),
+    ],
+)
+def test_slices_wrong_limit_rounding_method(limit_rounding_method, error_msg):
+    """Call overlap_slices with an invalid rounding method."""
+    with pytest.raises(ValueError, match=error_msg):
+        overlap_slices((5,), (3,), (0,), limit_rounding_method=limit_rounding_method)
+
+
 def test_extract_array_even_shape_rounding():
     """
     Test overlap_slices (via extract_array) for rounding with an
@@ -554,6 +632,25 @@ class TestCutout2D:
         c2.data[xy] = value
         assert data[yx] != value
 
+    @pytest.mark.parametrize(
+        "position, limit_rounding_method, expected_slices_original",
+        [
+            ((2, 2), np.ceil, (slice(1, 4), slice(1, 4))),
+            ((2, 2), np.floor, (slice(0, 3), slice(0, 3))),
+            ((1.9, 2.9), np.round, (slice(1, 4), slice(0, 3))),
+            ((2, 2), np.trunc, (slice(0, 3), slice(0, 3))),
+        ],
+    )
+    def test_limit_rounding_method(
+        self, position, limit_rounding_method, expected_slices_original
+    ):
+        c = Cutout2D(
+            self.data, position, (3, 3), limit_rounding_method=limit_rounding_method
+        )
+        assert c.data.shape == (3, 3)
+        assert c.slices_original == expected_slices_original
+        assert c.slices_cutout == (slice(0, 3), slice(0, 3))
+
     def test_to_from_large(self):
         position = (2, 2)
         c = Cutout2D(self.data, position, (3, 3))