diff --git a/cxx/isce3/core/LUT2d.cpp b/cxx/isce3/core/LUT2d.cpp
index 5da3a24f3..d24dfe141 100644
--- a/cxx/isce3/core/LUT2d.cpp
+++ b/cxx/isce3/core/LUT2d.cpp
@@ -164,6 +164,23 @@ eval(double y, const Eigen::Ref<const Eigen::VectorXd>& x) const
 {
     const auto n = x.size();
     Eigen::Matrix<T, Eigen::Dynamic, 1> out(n);
+
+    // Check bounds before parallel region to avoid exceptions in OpenMP threads
+    if (_boundsError && _haveData) {
+        for (long i = 0; i < n; ++i) {
+            if (!contains(y, x(i))) {
+                pyre::journal::error_t errorChannel("isce.core.LUT2d");
+                errorChannel
+                    << "Out of bounds LUT2d evaluation at " << y << " " << x(i)
+                    << pyre::journal::newline
+                    << " - bounds are " << _ystart << " "
+                    << _ystart + _dy * (_data.length() - 1.0) << " "
+                    << _xstart << " " << _xstart + _dx * (_data.width() - 1.0)
+                    << pyre::journal::endl;
+            }
+        }
+    }
+
     _Pragma("omp parallel for")
     for (long i = 0; i < n; ++i) {
         out(i) = eval(y, x(i));
diff --git a/tests/python/extensions/pybind/core/LUT2d.py b/tests/python/extensions/pybind/core/LUT2d.py
index 07ec08291..5ea37c24d 100644
--- a/tests/python/extensions/pybind/core/LUT2d.py
+++ b/tests/python/extensions/pybind/core/LUT2d.py
@@ -1,15 +1,17 @@
 #!/usr/bin/env python3
 
+import iscetest
+import journal
 import numpy as np
 
 import isce3.ext.isce3 as isce
-import iscetest
+
 
 def test_LUT2d():
     # Create LUT2d obj
     xvec = yvec = np.arange(-5.01, 5.01, 0.25)
     xx, yy = np.meshgrid(xvec, xvec)
-    M = np.sin(xx*xx + yy*yy)
+    M = np.sin(xx * xx + yy * yy)
     method = isce.core.DataInterpMethod.BIQUINTIC
     lut2d = isce.core.LUT2d(xvec, yvec, M, "biquintic")
     assert lut2d.interp_method == method
@@ -25,17 +27,17 @@ def test_LUT2d():
     assert lut2d.y_end == yvec[-1]
 
     # Load reference data
-    f_ref = iscetest.data + 'interpolator/data.txt'
+    f_ref = iscetest.data + "interpolator/data.txt"
     d_refs = np.loadtxt(f_ref)
-    
+
     # Loop over test points and check for error
     error = 0
     for d_ref in d_refs:
         z_test = lut2d.eval(d_ref[0], d_ref[1])
-        error += (d_ref[5] - z_test)**2
+        error += (d_ref[5] - z_test) ** 2
 
     n_pts = d_refs.shape[0]
-    assert error/n_pts < 0.058, f'pybind LUT2d failed: {error} > 0.058'
+    assert error / n_pts < 0.058, f"pybind LUT2d failed: {error} > 0.058"
 
     # check that we can set ref_value
     lut = isce.core.LUT2d()
@@ -44,3 +46,44 @@ def test_LUT2d():
     assert lut.ref_value == 1.0
     lut = isce.core.LUT2d(2.0)
     assert lut.ref_value == 2.0
+
+
+def test_bounds_error():
+    """Test that out-of-bounds evaluation raises exceptions.
+
+    Regression test for bug where vectorized .eval crashes Python instead of
+    raising a catchable exception when bounds_error=True.
+    """
+    import pytest
+
+    # Create a simple LUT with known bounds: 10 points in x, 2 points in y
+    x = np.linspace(0, 5, 10)
+    y = np.linspace(10, 20, 2)
+    z = np.vstack((np.linspace(100, 200, 10), np.linspace(100, 200, 10)))
+    lut2d = isce.core.LUT2d(x, y, z)
+    # Default is bounds_error=True
+    assert lut2d.bounds_error is True
+
+    y = 15.0
+    x = 2.0
+    lut2d.eval(y, x)  # y=1.0 is out of bounds (valid: 10-20)
+
+    # Test scalar out-of-bounds
+    x_oob = 200.0
+    with pytest.raises(journal.ApplicationError):
+        lut2d.eval(y, x_oob)
+
+    # Test vectorized out-of-bounds (should raise ApplicationError, not crash)
+    with pytest.raises(journal.ApplicationError):
+        lut2d.eval(y, np.array([x_oob, x_oob]))
+
+    # Test vectorized, all in-bounds
+    result = lut2d.eval(y, np.array([1.0, 2.0, 3.0]))
+    assert result.shape == (3,)
+
+    # Test with bounds_error=False to avoid raising exceptions
+    lut2d.bounds_error = False
+    result = lut2d.eval(y, x_oob)  # Should clamp and return value
+    assert result == 200.0
+    result = lut2d.eval(y, np.array([x_oob, x_oob]))  # Should clamp and return values
+    assert np.allclose(result, np.array([200, 200]))