diff --git a/gwcs/wcs.py b/gwcs/wcs.py
index e44bde42..baaf5983 100644
--- a/gwcs/wcs.py
+++ b/gwcs/wcs.py
@@ -551,11 +551,16 @@ def _call_backward(self, *args, with_bounding_box=True, fill_value=np.nan, **kwa
         return result
 
     def outside_footprint(self, world_arrays):
+        world_eps = 5.0 * np.finfo(np.float32).eps
+        eps_p = 1.0 + world_eps
+        eps_m = 1.0 - world_eps
+
         world_arrays = list(world_arrays)
 
         axes_types = set(self.output_frame.axes_type)
         axes_phys_types = self.world_axis_physical_types
         footprint = self.footprint()
+
         not_numerical = False
         if utils.is_high_level(world_arrays[0], low_level_wcs=self):
             not_numerical = True
@@ -579,9 +584,10 @@ def outside_footprint(self, world_arrays):
                             axis_range[axis_range < d].max(),
                             axis_range[axis_range > d].min()
                         ]
-                        outside = (coord >= range[0]) & (coord < range[1])
+                        outside = (coord > range[0] * eps_p) & (coord < range[1] * eps_m)
                 else:
-                    outside = (coord < range[0]) | (coord > range[1])
+                    outside = (coord < range[0] * eps_m) | (coord > range[1] * eps_p)
+
                 if np.any(outside):
                     if np.isscalar(coord):
                         coord = np.nan