fix(geospatial): fix VertexGrid layer search and boundary point handling

Fixes two critical issues in GeospatialIndex:

1. VertexGrid layer search: Was always returning layer 0. Now properly
   searches through all layers using zv[lay, icell2d] for top and
   zv[lay+1, icell2d] for bottom elevations.

2. Boundary point handling: Points exactly on shared vertices/edges
   could match multiple cells. Implemented max(matching_cells) tiebreaker
   to match searchsorted() behavior with side="right" for consistency.

Performance impact: Checks up to 5 candidates instead of returning first
match, adding minimal overhead while ensuring correctness for edge cases.

Fixes:
- test_vertex_xyz_intersect (layer 0 vs layer 2 mismatch)
- test_intersection (cell 267 vs 268 boundary point)

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

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

Author: adacovsk

flopy/utils/geospatial_index.py

@@ -388,26 +388,40 @@ 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 matching cells for boundary point tiebreaker
+            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)
+                        # For performance, stop after finding enough matches
+                        if len(matching_cells) >= 5:
+                            break
                 else:
                     # 2D: use ConvexHull test
                     if self._point_in_cell_vectorized(point[:2], cellid):
-                        # Found 2D cell, now check z if provided
+                        # For z-search, handle differently
                         if z is None:
-                            results[i] = cellid
-                            break
+                            matching_cells.append(cellid)
+                            # For performance, stop after finding enough matches
+                            if len(matching_cells) >= 5:
+                                break
                         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)
+                                # For performance, stop after finding enough matches
+                                if len(matching_cells) >= 5:
+                                    break
+
+            # If multiple cells match (boundary point), choose highest index
+            # This matches searchsorted() behavior with side="right"
+            if matching_cells:
+                results[i] = max(matching_cells)
 
         return results
 
@@ -554,20 +568,19 @@ def _find_layer_for_z(self, icell2d, z):
 
         # For VertexGrid: same 2D geometry for all layers
         if self.grid.grid_type == "vertex":
-            # Build array of cell indices for this 2D cell across all layers
-            # For VertexGrid, it's just the same icell2d index repeated
-            z_top = zv[0, icell2d]
-            z_bot = zv[1, icell2d]
-            z_min = min(z_top, z_bot)
-            z_max = max(z_top, z_bot)
-
-            # Check if z is within bounds
-            if z_min <= z <= z_max:
-                # Find which layer by proportional division
-                # For now, just return first layer (layer index)
-                # Note: VertexGrid may have varying layer thicknesses
-                # TODO: This assumes uniform z for all layers - may need refinement
-                return 0 if z_min <= z <= z_max else None
+            # Search through all layers to find which contains z
+            # For VertexGrid: zv[lay, icell2d] = top, zv[lay+1, icell2d] = bottom
+            for lay in range(self.grid.nlay):
+                z_top = zv[lay, icell2d]
+                z_bot = zv[lay + 1, icell2d]
+                z_min = min(z_top, z_bot)
+                z_max = max(z_top, z_bot)
+
+                if z_min <= z <= z_max:
+                    # Found the layer! Return the layer index
+                    return lay
+
+            # z not found in any layer
             return None
 
         # For UnstructuredGrid: compute 3D cell index using vectorized search