diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py
index f1e34cb9..cc0e3a14 100644
--- a/gwcs/tests/test_coordinate_systems.py
+++ b/gwcs/tests/test_coordinate_systems.py
@@ -112,6 +112,7 @@ def test_bare_baseframe():
             output_frame=frame,
             input_frame=cf.CoordinateFrame(1, "PIXEL", (0,), unit=(u.pix,), name="detector_frame")
             )
+    #w.bounding_box = (0, 9)
     assert u.allclose(w.world_to_pixel(0*u.km), 0)
 
 
diff --git a/gwcs/wcs.py b/gwcs/wcs.py
index 6e5a4a64..3d281f01 100644
--- a/gwcs/wcs.py
+++ b/gwcs/wcs.py
@@ -480,9 +480,13 @@ def invert(self, *args, **kwargs):
 
         with_bounding_box = kwargs.pop('with_bounding_box', True)
         fill_value = kwargs.pop('fill_value', np.nan)
+        akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
+
+        if with_bounding_box and self.bounding_box is not None:
+            result = self.outside_footprint(args)
 
         if btrans is not None:
-            akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
+            #akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS}
             result = btrans(*args, **akwargs)
         else:
             result = self.numerical_inverse(*args, **kwargs, with_units=with_units)
@@ -498,6 +502,26 @@ def invert(self, *args, **kwargs):
                 return self.input_frame.coordinates(*result)
         else:
             return result
+        
+    def outside_footprint(self, world_arrays):
+        for axis in world_arrays:
+            if np.isscalar(world_arrays) or self.output_frame.naxes == 1:
+                world_arrays = [world_arrays]
+        world_arrays = list(world_arrays)
+        footprint = self.footprint()
+        for idim, coord in enumerate(world_arrays):
+            axis_range = footprint[:, idim]
+            range = [axis_range.min(), axis_range.max()]
+            outside = (coord < range[0]) | (coord > range[1])
+            if np.any(outside):
+                if np.isscalar(coord):
+                    coord = np.nan
+                else:
+                    coord[outside] = np.nan
+                world_arrays[idim] = coord
+
+        return world_arrays
+
 
     def out_of_bounds(self, pixel_arrays, fill_value=np.nan):
         if np.isscalar(pixel_arrays) or self.input_frame.naxes == 1:
@@ -1401,6 +1425,11 @@ def footprint(self, bounding_box=None, center=False, axis_type="all"):
 
         """
         def _order_clockwise(v):
+            # if self.input_frame.naxes == 1:
+            #     bb = self.bounding_box.bounding_box()
+            #     if isinstance(bb[0], u.Quantity):
+            #         bb = [v.value for v in bb] * bb[0].unit
+            #     return (bb,)
             return np.asarray([[v[0][0], v[1][0]], [v[0][0], v[1][1]],
                                [v[0][1], v[1][1]], [v[0][1], v[1][0]]]).T
 
@@ -1418,6 +1447,7 @@ def _order_clockwise(v):
         else:
             vertices = np.array(list(itertools.product(*bb))).T
 
+        # workaround an issue with bbox with quantity, interval needs to be a cquantity, not a list of quantities
         if center:
             vertices = utils._toindex(vertices)