Possible fix for world-to-pixel when out of bounding box

CHANGES.rst

@@ -1,6 +1,8 @@
 0.22.0 (unreleased)
 -------------------
 
+- Return nan for world to pixel transformations outside of the
+  bounding box. [#497]
 
 0.21.0 (2024-03-10)
 -------------------

gwcs/api.py

@@ -135,8 +135,10 @@ def world_to_pixel_values(self, *world_arrays):
         world_arrays = self._add_units_input(world_arrays, self.backward_transform, self.output_frame)
 
         result = self.invert(*world_arrays, with_units=False)
+        result = self._remove_quantity_output(result, self.input_frame)
+
+        return self._out_of_bounding_box_to_nan(result, world_arrays)
 
-        return self._remove_quantity_output(result, self.input_frame)
 
     def world_to_array_index_values(self, *world_arrays):
         """
@@ -322,6 +324,8 @@ def world_to_pixel(self, *world_objects):
             if not utils.isnumerical(result):
                 result = result.value
 
+        result = self._out_of_bounding_box_to_nan(result, world_objects)
+
         return result
 
     def world_to_array_index(self, *world_objects):

gwcs/wcs.py

@@ -396,7 +396,11 @@ def in_image(self, *args, **kwargs):
         kwargs['fill_value'] = np.nan
 
         coords = self.invert(*args, **kwargs)
+        result = self._pixel_coords_in_image(coords, args)
 
+        return result
+
+    def _pixel_coords_in_image(self, coords, args):
         result = np.isfinite(coords)
         if self.input_frame.naxes > 1:
             result = np.all(result, axis=0)
@@ -419,6 +423,20 @@ def in_image(self, *args, **kwargs):
 
             elif result:
                 result = all([(c >= x1) and (c <= x2) for c, (x1, x2) in zip(coords, self.bounding_box)])
+        return result
+
+    def _out_of_bounding_box_to_nan(self, coords, args):
+        # world coordinates outside the bounding_box are set to nan:
+        within_bounding_box = self._pixel_coords_in_image(coords, args)
+        result = np.where(
+            within_bounding_box,
+            coords,
+            np.nan
+        )
+
+        if isinstance(result, np.ndarray) and result.shape == ():
+            # scalar results should be floats:
+            result = float(result)
 
         return result