Fix ravel model behaviour to properly clip OOB coordinates

changelog/751.bugfix.rst

@@ -0,0 +1 @@
+Fix coordinate clipping behaviour of Ravel model.

dkist/dataset/tests/test_dataset.py

@@ -14,6 +14,7 @@
 from dkist.dataset import Dataset, TiledDataset, load_dataset
 from dkist.io import DKISTFileManager
 from dkist.utils.exceptions import DKISTDeprecationWarning
+from dkist.wcs.models import Ravel
 
 
 @pytest.fixture
@@ -213,3 +214,28 @@ def test_file_slicing_without_dummy_axis(dataset_5d):
     assert len(ds[0, 0].files) == np.prod(shape[2])
     assert len(ds[0, 0, 0].files) == 1
     assert len(ds[0, 0, 0, 0].files) == 1
+
+
+@pytest.mark.parametrize("ndim", [pytest.param(2, id="2D"), pytest.param(3, id="3D")])
+@pytest.mark.parametrize("has_units", [pytest.param(True, id="With Units"), pytest.param(False, id="Without Units")])
+@pytest.mark.parametrize("input_type", [pytest.param("array", id="Array Inputs"), pytest.param("scalar", id="Scalar Inputs")])
+def test_ravel_oob_nan(ndim, has_units, input_type):
+    array_shape = tuple(5 for _ in range(ndim))
+    ravel = Ravel(array_shape, order="C")
+    units = u.pix
+    if input_type == "array":
+        # One in-bounds, two out-of-bounds (below and above)
+        coords = [np.array([2, -1, 10])] + [np.array([2, 2, 2])] * (ndim - 1)
+    else:
+        coords = [(-1)] + [2] * (ndim - 1)
+    if has_units:
+        coords = coords * units
+    result = ravel(*coords)
+    # First element (in-bounds) must be finite, OOB elements must be NaN
+    if input_type == "array":
+        out = result.value if has_units else result
+        assert np.isfinite(out[0])
+        assert np.isnan(out[1])
+        assert np.isnan(out[2])
+    else:
+        assert np.isnan(result)

dkist/wcs/models.py

@@ -739,17 +739,19 @@ def evaluate(self, *inputs_):
         # round the index values, but clip them if they exceed the array bounds
         # the bounds are one less than the shape dimension value
         array_bounds = np.array(self.array_shape) - 1
-        rounded_inputs = np.clip(np.rint(input_values).astype(int), None, array_bounds[:, np.newaxis])
-        result = np.ravel_multi_index(rounded_inputs, self.array_shape, order=self.order).astype(float)
+        rounded_inputs = np.rint(input_values).astype(int)
+        result = np.ravel_multi_index(rounded_inputs, self.array_shape, order=self.order, mode="clip").astype(float)
         index = 0 if self.order == "F" else -1
         # Adjust the result to allow a fractional part for interpolation in Tabular1D
         fraction = input_values[index] - rounded_inputs[index]
         result += fraction
+        oob = np.logical_or((rounded_inputs < 0), (rounded_inputs > array_bounds[:, np.newaxis])).any(axis=0)
+        result[oob] = np.nan
         # Put the units back if they were there...
         if has_units:
             result = result * u.pix
         else:
-            result = np.array([result])
+            result = np.array(result)
         return result
 
     @property