diff --git a/changelog/751.bugfix.rst b/changelog/751.bugfix.rst new file mode 100644 index 00000000..5231faa3 --- /dev/null +++ b/changelog/751.bugfix.rst @@ -0,0 +1 @@ +Fix coordinate clipping behaviour of Ravel model. diff --git a/dkist/dataset/tests/test_dataset.py b/dkist/dataset/tests/test_dataset.py index 678567bb..40394e36 100644 --- a/dkist/dataset/tests/test_dataset.py +++ b/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) diff --git a/dkist/wcs/models.py b/dkist/wcs/models.py index 5440ab51..f19fb75b 100755 --- a/dkist/wcs/models.py +++ b/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