diff --git a/gwcs/wcs.py b/gwcs/wcs.py
index db0b1a80..a52dc06e 100644
--- a/gwcs/wcs.py
+++ b/gwcs/wcs.py
@@ -15,7 +15,7 @@
                                      Sky2Pix_TAN)
 from astropy.modeling.parameters import _tofloat
 from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales
-from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values
+from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values, values_to_high_level_objects
 
 from astropy import units as u
 from scipy import linalg, optimize
@@ -494,7 +494,9 @@ def outside_footprint(self, world_arrays):
 
         axes_types = set(self.output_frame.axes_type)
         footprint = self.footprint()
+        not_numerical = False
         if not utils.isnumerical(world_arrays[0]):
+            not_numerical = True
             world_arrays = high_level_objects_to_values(*world_arrays, low_level_wcs=self)
         for axtyp in axes_types:
             ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp))
@@ -513,7 +515,8 @@ def outside_footprint(self, world_arrays):
                     else:
                         coord[outside] = np.nan
                     world_arrays[idim] = coord
-
+        if not_numerical:
+            world_arrays = values_to_high_level_objects(*world_arrays, low_level_wcs=self)
         return world_arrays