Return NaN for out of bound conversions in FITS WCS and update APE 14…

… to say that implementations should return NaN rather than may.

astropy/wcs/wcsapi/fitswcs.py

@@ -332,7 +332,26 @@ def axis_correlation_matrix(self):
 
         return matrix
 
+    def _out_of_bounds_to_nan(self, pixel_arrays):
+        if self.pixel_bounds is not None:
+            pixel_arrays = list(pixel_arrays)
+            for idim in range(self.pixel_n_dim):
+                if self.pixel_bounds[idim] is None:
+                    continue
+                out_of_bounds = (pixel_arrays[idim] < self.pixel_bounds[idim][0]) | (
+                    pixel_arrays[idim] > self.pixel_bounds[idim][1]
+                )
+                if np.any(out_of_bounds):
+                    pix = pixel_arrays[idim].astype(float, copy=True)
+                    if np.isscalar(pix):
+                        pix = np.nan
+                    else:
+                        pix[out_of_bounds] = np.nan
+                    pixel_arrays[idim] = pix
+        return pixel_arrays
+
     def pixel_to_world_values(self, *pixel_arrays):
+        pixel_arrays = self._out_of_bounds_to_nan(pixel_arrays)
         world = self.all_pix2world(*pixel_arrays, 0)
         return world[0] if self.world_n_dim == 1 else tuple(world)
 
@@ -350,6 +369,8 @@ def world_to_pixel_values(self, *world_arrays):
                 lambda *args: e.best_solution, "input", *world_arrays, 0
             )
 
+        pixel = self._out_of_bounds_to_nan(pixel)
+
         return pixel[0] if self.pixel_n_dim == 1 else tuple(pixel)
 
     @property

astropy/wcs/wcsapi/low_level_api.py

@@ -64,7 +64,7 @@ def pixel_to_world_values(self, *pixel_arrays):
         `~astropy.wcs.wcsapi.BaseLowLevelWCS.world_n_dim` scalars or arrays in units given by
         `~astropy.wcs.wcsapi.BaseLowLevelWCS.world_axis_units`. Note that pixel coordinates are
         assumed to be 0 at the center of the first pixel in each dimension. If a
-        pixel is in a region where the WCS is not defined, NaN can be returned.
+        pixel is in a region where the WCS is not defined, NaN should be returned.
         The coordinates should be specified in the ``(x, y)`` order, where for
         an image, ``x`` is the horizontal coordinate and ``y`` is the vertical
         coordinate.
@@ -99,7 +99,7 @@ def world_to_pixel_values(self, *world_arrays):
         `~astropy.wcs.wcsapi.BaseLowLevelWCS.pixel_n_dim` scalars or arrays. Note that pixel
         coordinates are assumed to be 0 at the center of the first pixel in each
         dimension. If a world coordinate does not have a matching pixel
-        coordinate, NaN can be returned.  The coordinates should be returned in
+        coordinate, NaN should be returned.  The coordinates should be returned in
         the ``(x, y)`` order, where for an image, ``x`` is the horizontal
         coordinate and ``y`` is the vertical coordinate.
 
@@ -274,6 +274,9 @@ def pixel_bounds(self):
         within a certain range of pixel values, for example when defining a
         WCS that includes fitted distortions. This is an optional property,
         and it should return `None` if a shape is not known or relevant.
+
+        The bounds can be a mix of values along dimensions where bounds exist,
+        and None for other dimensions, e.g. ``[(xmin, xmax), None]``.
         """
         return None
 

astropy/wcs/wcsapi/tests/test_fitswcs.py

@@ -1417,3 +1417,87 @@ def test_pixel_to_world_stokes(wcs_polarized):
     assert isinstance(world[2], StokesCoord)
     assert_array_equal(world[2], [1, 2, 3, 4])
     assert_array_equal(world[2].symbol, ["I", "Q", "U", "V"])
+
+
+@pytest.mark.parametrize("direction", ("world_to_pixel", "pixel_to_world"))
+def test_out_of_bounds(direction):
+    # Make sure that we correctly deal with any out-of-bound values in the
+    # low-level API.
+
+    wcs = WCS(naxis=2)
+    wcs.wcs.crpix = (1, 1)
+    wcs.wcs.set()
+
+    func = (
+        wcs.world_to_pixel_values
+        if direction == "world_to_pixel"
+        else wcs.pixel_to_world_values
+    )
+
+    xp = np.arange(5) + 1
+    yp = np.arange(5) + 1
+
+    # Before setting bounds
+
+    # Scalars
+    xw, yw = func(xp[0], yp[0])
+    assert_array_equal(xw, 1)
+    assert_array_equal(yw, 1)
+
+    # Arrays
+    xw, yw = func(xp, yp)
+    assert_array_equal(xw, [1, 2, 3, 4, 5])
+    assert_array_equal(yw, [1, 2, 3, 4, 5])
+
+    # Mixed
+    xw, yw = func(xp[0], yp)
+    assert_array_equal(xw, 1)
+    assert_array_equal(yw, [1, 2, 3, 4, 5])
+
+    # Setting bounds on one dimension
+
+    wcs.pixel_bounds = [(-0.5, 3.5), None]
+
+    # Scalars
+
+    xw, yw = func(xp[0], yp[0])
+    assert_array_equal(xw, 1)
+    assert_array_equal(yw, 1)
+
+    xw, yw = func(xp[-1], yp[-1])
+    assert_array_equal(xw, np.nan)
+    assert_array_equal(yw, 5)
+
+    # Arrays
+    xw, yw = func(xp, yp)
+    assert_array_equal(xw, [1, 2, 3, np.nan, np.nan])
+    assert_array_equal(yw, [1, 2, 3, 4, 5])
+
+    # Mixed
+    xw, yw = func(xp[-1], yp)
+    assert_array_equal(xw, np.nan)
+    assert_array_equal(yw, [1, 2, 3, 4, 5])
+
+    # Setting bounds on both dimensions
+
+    wcs.pixel_bounds = [(-0.5, 3.5), (2.5, 5.5)]
+
+    # Scalars
+
+    xw, yw = func(xp[0], yp[0])
+    assert_array_equal(xw, 1)
+    assert_array_equal(yw, np.nan)
+
+    xw, yw = func(xp[-1], yp[-1])
+    assert_array_equal(xw, np.nan)
+    assert_array_equal(yw, 5)
+
+    # Arrays
+    xw, yw = func(xp, yp)
+    assert_array_equal(xw, [1, 2, 3, np.nan, np.nan])
+    assert_array_equal(yw, [np.nan, np.nan, 3, 4, 5])
+
+    # Mixed
+    xw, yw = func(xp[-1], yp)
+    assert_array_equal(xw, np.nan)
+    assert_array_equal(yw, [np.nan, np.nan, 3, 4, 5])