fix(geospatial): add epsilon tolerance for boundary points

Add configurable epsilon parameter to GeospatialIndex to handle points
on or near cell boundaries. The epsilon is used for:
- Expanding bounding box tests (bbox ± epsilon)
- ConvexHull hyperplane distance tests (distance <= epsilon)

Default epsilon=1e-6 works well for grids with coordinates up to ~10,000 units.

This fixes errors where valid points near boundaries were incorrectly
rejected as outside the model area, such as:
- point (0.1, 8999.9, -16.67) in DISU grids
- point (0.1, 10499.9, 310.0) in larger grids

The old matplotlib Path.contains_point() approach used radius=1e-9 to
include edge points. The new ConvexHull approach uses epsilon=1e-6 which
provides similar edge inclusion while being more robust for larger grids.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude <[email protected]>

Author: adacovsk

flopy/utils/geospatial_index.py

@@ -26,11 +26,17 @@ class GeospatialIndex:
     ----------
     grid : Grid
         Any FloPy grid object (StructuredGrid, VertexGrid, UnstructuredGrid)
+    epsilon : float, optional
+        Tolerance for point-in-cell tests. Used for both bounding box
+        expansion and ConvexHull hyperplane distance tests. Default is 1e-6,
+        which works well for grids with coordinates up to ~10,000 units.
 
     Attributes
     ----------
     grid : Grid
         The grid object this index was built for
+    epsilon : float
+        Tolerance for geometric tests
     tree : scipy.spatial.KDTree
         KD-tree of cell centroids + vertices for fast spatial queries
     point_to_cell : ndarray
@@ -55,16 +61,20 @@ class GeospatialIndex:
     >>> cellids = index.query_points(x=[1.5, 5.5, 9.5], y=[1.5, 5.5, 9.5])
     """
 
-    def __init__(self, grid):
+    def __init__(self, grid, epsilon=1e-6):
         """
         Build geospatial index for a model grid.
 
         Parameters
         ----------
         grid : Grid
             Any FloPy grid (StructuredGrid, VertexGrid, UnstructuredGrid)
+        epsilon : float, optional
+            Tolerance for point-in-cell tests (default 1e-6).
+            Used for bounding box expansion and ConvexHull tests.
         """
         self.grid = grid
+        self.epsilon = epsilon
         self._build_index()
 
     def _build_index(self):
@@ -357,17 +367,20 @@ def _point_in_cell_vectorized(self, point, cellid):
         bool
             True if point is inside cell
         """
-        # Fast bounding box rejection
+        # Fast bounding box rejection with epsilon tolerance for edge cases
+        # Expand bounding box slightly to handle points on boundaries
         bbox = self.bounding_boxes[cellid]
-        if not (bbox[0] <= point[0] <= bbox[1] and bbox[2] <= point[1] <= bbox[3]):
+        if not (bbox[0] - self.epsilon <= point[0] <= bbox[1] + self.epsilon and
+                bbox[2] - self.epsilon <= point[1] <= bbox[3] + self.epsilon):
             return False
 
         # Precise ConvexHull test
         hull_eq = self.hull_equations[cellid]
         if hull_eq is not None:
-            # Vectorized hyperplane test: Ax + By + C <= tolerance
+            # Vectorized hyperplane test: Ax + By + C <= epsilon
+            # Tolerance allows points on edges to be included
             distances = point @ hull_eq[:, :-1].T + hull_eq[:, -1]
-            return np.all(distances <= 1e-6)
+            return np.all(distances <= self.epsilon)
         else:
             # Degenerate geometry - should rarely happen
             return False