diff --git a/CHANGES.rst b/CHANGES.rst
index c8f89457bc3..e2b625443d0 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -137,6 +137,10 @@ resample_spec
 
 - Populate the wavelength array in resampled `Slit` and `MultiSlit` models. [#8374]
 
+- Fix a bug resulting in large WCS errors in the resampled image's WCS
+  under certain circumstances (depending on the orientation of the slit on the
+  sky). [#8511]
+
 residual_fringe
 ---------------
 
diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py
index fcd49a93029..71413ec9a3a 100644
--- a/jwst/resample/resample_spec.py
+++ b/jwst/resample/resample_spec.py
@@ -408,7 +408,7 @@ def build_interpolated_output_wcs(self, refmodel=None):
                 undist2sky1 = tan | native2celestial
                 # Filter out RuntimeWarnings due to computed NaNs in the WCS
                 with warnings.catch_warnings():
-                    warnings.simplefilter("ignore", RuntimeWarning) #was ignore. need to make more specific
+                    warnings.simplefilter("ignore", RuntimeWarning)  # was ignore. need to make more specific
                 # at this center of slit find x,y tangent projection - x_tan, y_tan
                 x_tan, y_tan = undist2sky1.inverse(ra, dec)
 
@@ -498,7 +498,7 @@ def build_interpolated_output_wcs(self, refmodel=None):
         # mapping.inverse uses pix_to_ytan.inverse.  The auto definition
         # of mapping.inverse is to use the 2nd spatial coordinate, i.e. Dec.
 
-        swap_xy = np.isclose(pix_to_xtan.slope, 0, atol=1e-8)
+        swap_xy = abs(pix_to_xtan.slope) < abs(pix_to_ytan.slope)
         if swap_xy:
             # Account for vertical or horizontal dispersion on detector
             mapping.inverse = Mapping((2, 1) if spatial_axis else (1, 2))