feat(geospatial): add tie-breaking for boundary points

Implements grid-specific tie-breaking rules when points fall on cell
boundaries and match multiple cells. This ensures deterministic cell
selection and prevents particle tracking issues.

Changes:
- GeospatialIndex: Collect all matching cells and apply tie-breaker
- StructuredGrid: Use searchsorted with side='left'/'right' for tie-breaking
- Add boundary detection methods for 2D and 3D points
- Implement _apply_tiebreaker() with grid-specific rules:
  - StructuredGrid: lowest row, then lowest column
  - VertexGrid/UnstructuredGrid: lowest cell ID

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

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

Author: adacovsk

flopy/discretization/structuredgrid.py

@@ -1029,17 +1029,19 @@ def intersect(self, x, y, z=None, local=False, forgive=False):
         cols = np.full(n_points, np.nan, dtype=float)
 
         # Vectorized column finding using searchsorted
-        # For each x[i], find rightmost edge where x[i] > edge
-        cols_valid = np.searchsorted(xe, x, side="right") - 1
+        # When point is on edge, use lower column index (tie-breaking rule)
+        # side="left" ensures x==edge goes to the cell on the left
+        cols_valid = np.searchsorted(xe, x, side="left") - 1
         # searchsorted returns index in [0, len(xe)], but valid cols are [0, ncol-1]
         cols_mask = (cols_valid >= 0) & (cols_valid < self.ncol)
         cols[cols_mask] = cols_valid[cols_mask]
 
         # Vectorized row finding using searchsorted
-        # For each y[i], find rightmost edge where y[i] < edge
+        # When point is on edge, use lower row index (tie-breaking rule)
         # Since ye is in descending order (top to bottom), we need to flip it
+        # side="right" on flipped array ensures y==edge goes to the cell below
         ye_flipped = ye[::-1]
-        rows_flipped = np.searchsorted(ye_flipped, y, side="left")
+        rows_flipped = np.searchsorted(ye_flipped, y, side="right")
         rows_valid = len(ye) - 1 - rows_flipped
         rows_mask = (rows_valid >= 0) & (rows_valid < self.nrow)
         rows[rows_mask] = rows_valid[rows_mask]

flopy/utils/geospatial_index.py

@@ -391,26 +391,30 @@ def query_points(self, x, y, z=None, k=None):
             # Query KD-tree adaptively to get k unique cells
             candidates = self._get_k_unique_cells(point, k)
 
+            # Collect all matching cells for tie-breaking
+            matching_cells = []
+
             # Test candidates in order (nearest first)
             for cellid in candidates:
                 if self.is_3d:
                     # 3D: check 3D bounding box
                     if self._point_in_cell_3d(point, cellid):
-                        results[i] = cellid
-                        break
+                        matching_cells.append(cellid)
                 else:
                     # 2D: use ConvexHull test
                     if self._point_in_cell_vectorized(point[:2], cellid):
                         # Found 2D cell, now check z if provided
                         if z is None:
-                            results[i] = cellid
-                            break
+                            matching_cells.append(cellid)
                         else:
                             # Search through layers to find the right z
                             cell_3d = self._find_layer_for_z(cellid, z[i])
                             if cell_3d is not None:
-                                results[i] = cell_3d
-                                break
+                                matching_cells.append(cell_3d)
+
+            # Apply grid-specific tie-breaking when multiple cells match
+            if len(matching_cells) > 0:
+                results[i] = self._apply_tiebreaker(matching_cells)
 
         return results
 
@@ -534,6 +538,87 @@ def _point_in_cell_3d(self, point, cellid):
             and bbox[4] - self.epsilon <= point[2] <= bbox[5] + self.epsilon
         )
 
+    def _is_on_boundary_2d(self, point, cellid):
+        """
+        Check if a 2D point is on the boundary of a cell.
+
+        A point is considered on the boundary if it's very close (< 1e-9) to any edge.
+        Uses a much tighter tolerance than epsilon to detect true boundary cases.
+        """
+        bbox = self.bounding_boxes[cellid]
+        boundary_tol = 1e-9  # Much tighter than epsilon for exact boundary detection
+        # Check if point is on any bounding box edge
+        on_x_edge = (
+            abs(point[0] - bbox[0]) < boundary_tol or
+            abs(point[0] - bbox[1]) < boundary_tol
+        )
+        on_y_edge = (
+            abs(point[1] - bbox[2]) < boundary_tol or
+            abs(point[1] - bbox[3]) < boundary_tol
+        )
+        return on_x_edge or on_y_edge
+
+    def _is_on_boundary_3d(self, point, cellid):
+        """
+        Check if a 3D point is on the boundary of a cell.
+
+        A point is considered on the boundary if it's very close (< 1e-9) to any face.
+        Uses a much tighter tolerance than epsilon to detect true boundary cases.
+        """
+        bbox = self.bounding_boxes_3d[cellid]
+        boundary_tol = 1e-9  # Much tighter than epsilon for exact boundary detection
+        # Check if point is on any bounding box face
+        on_x_face = (
+            abs(point[0] - bbox[0]) < boundary_tol or
+            abs(point[0] - bbox[1]) < boundary_tol
+        )
+        on_y_face = (
+            abs(point[1] - bbox[2]) < boundary_tol or
+            abs(point[1] - bbox[3]) < boundary_tol
+        )
+        on_z_face = (
+            abs(point[2] - bbox[4]) < boundary_tol or
+            abs(point[2] - bbox[5]) < boundary_tol
+        )
+        return on_x_face or on_y_face or on_z_face
+
+    def _apply_tiebreaker(self, matching_cells):
+        """
+        Apply grid-specific tie-breaking when multiple cells match.
+
+        Implements the same tie-breaking logic as the original grid implementations:
+        - StructuredGrid: cell with lowest row, then lowest column
+        - VertexGrid/UnstructuredGrid: cell with lowest cell ID
+
+        Parameters
+        ----------
+        matching_cells : list
+            List of cell indices that all contain the point
+
+        Returns
+        -------
+        cellid : int
+            The selected cell after tie-breaking
+        """
+        if len(matching_cells) == 1:
+            return matching_cells[0]
+
+        if self.grid.grid_type == "structured":
+            # For structured grids: choose cell with lowest row, then lowest column
+            # Convert cell indices to (row, col) tuples for comparison
+            rows_cols = []
+            for cellid in matching_cells:
+                row = cellid // self.grid.ncol
+                col = cellid % self.grid.ncol
+                rows_cols.append((row, col, cellid))
+
+            # Sort by row first, then column
+            rows_cols.sort(key=lambda x: (x[0], x[1]))
+            return rows_cols[0][2]  # Return the cellid
+        else:
+            # For vertex/unstructured grids: choose cell with lowest cell ID
+            return min(matching_cells)
+
     def _find_layer_for_z(self, icell2d, z):
         """
         Find 3D cell index for a 2D cell and z-coordinate.