feat(geospatial): add KD-tree spatial indexing with 3D support

Implement GeospatialIndex class using KD-tree and ConvexHull for fast
spatial queries on unstructured grids. Provides massive performance
improvements for grid.intersect() operations on VertexGrid and
UnstructuredGrid with up to 414x speedup for large grids.

**Architecture:**
- 2D mode: (x,y) KD-tree + ConvexHull for point-in-polygon tests
- 3D mode: (x,y,z) KD-tree + 3D bounding boxes for grid_varies_by_layer
- Automatic mode detection based on grid properties
- Epsilon tolerance (default 1e-6) for boundary point handling
- Vectorized z-coordinate search using numpy masks

**Performance Benchmarks:**
- VertexGrid: 1.2x-348x speedup (best for 1000+ cells)
- UnstructuredGrid 2D: 1.1x-352x speedup (best for 1000+ cells)
- UnstructuredGrid 3D: 1.1x-414x speedup (best for 2000+ cells)
- StructuredGrid: 0.10x-0.58x (slower, keeps using searchsorted)

**Key Features:**
- Pre-computed ConvexHull equations for vectorized containment tests
- Bounding box rejection for fast filtering
- KD-tree nearest neighbor search to find candidate cells
- Supports both 2D and 3D spatial indexing
- Handles grid_varies_by_layer for true 3D grids
- Epsilon tolerance for robust boundary point handling

**Implementation Details:**
- GeospatialIndex class in flopy/utils/geospatial_index.py
- _build_index() for 2D (x,y) point indexing with centroids + vertices
- _build_3d_index() for 3D (x,y,z) point indexing
- _precompute_hulls() for 2D bounding boxes and ConvexHull equations
- _precompute_3d_bounds() for 3D bounding boxes
- _point_in_cell_vectorized() for 2D containment tests
- _point_in_cell_3d() for 3D bounding box containment
- _find_layer_for_z() vectorized layer search (no loops)
- query_point() and query_points() public API

**Grid Integration:**
- StructuredGrid: continues using searchsorted() (faster)
- VertexGrid: uses GeospatialIndex for 2D + z-search
- UnstructuredGrid: uses GeospatialIndex for both 2D and 3D modes
- Automatic fallback to old methods for small grids

**Testing:**
- autotest/test_edge_cases.py: boundary conditions and edge cases
- autotest/test_geospatial_benchmark.py: performance validation

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

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

Author: adacovsk

autotest/test_edge_cases.py

@@ -0,0 +1,213 @@
+"""
+Edge case tests for GeospatialIndex to ensure robustness.
+"""
+
+import numpy as np
+import pytest
+from scipy.spatial import Delaunay
+
+from flopy.discretization import VertexGrid
+from flopy.utils.geospatial_index import GeospatialIndex
+
+
+def test_thin_sliver_cell():
+    """
+    Test that GeospatialIndex can find points in very thin "sliver" cells
+    where the centroid might be far from the actual cell location.
+    """
+    # Create a grid with a very thin sliver cell
+    # This tests the centroid+vertices KD-tree approach
+    np.random.seed(42)
+
+    # Create base random points
+    n_points = 15
+    x_verts = np.random.uniform(0, 100, n_points).tolist()
+    y_verts = np.random.uniform(0, 100, n_points).tolist()
+
+    # Add vertices for a very thin vertical sliver at x=50
+    sliver_indices = []
+    for i in range(4):
+        idx = len(x_verts)
+        sliver_indices.append(idx)
+        x_verts.append(50.0 + i * 0.05)  # Very thin: 0.15 units wide
+        y_verts.append(i * 33.33)  # Tall: 100 units high
+
+    # Create Delaunay triangulation
+    points = np.column_stack([x_verts, y_verts])
+    tri = Delaunay(points)
+
+    # Build VertexGrid
+    vertices = [[i, x_verts[i], y_verts[i]] for i in range(len(x_verts))]
+    cell2d = []
+    for i, simplex in enumerate(tri.simplices):
+        # Calculate centroid
+        cell_x = np.mean([x_verts[j] for j in simplex])
+        cell_y = np.mean([y_verts[j] for j in simplex])
+        cell2d.append([i, cell_x, cell_y, len(simplex)] + list(simplex))
+
+    ncells = len(cell2d)
+    grid = VertexGrid(
+        vertices=vertices,
+        cell2d=cell2d,
+        top=np.ones(ncells) * 10.0,
+        botm=np.zeros(ncells),
+    )
+
+    # Build index
+    index = GeospatialIndex(grid)
+
+    # Test points in/near the sliver region
+    test_points = [
+        (50.025, 50.0),  # Should be in a sliver cell
+        (50.075, 25.0),  # Should be in a sliver cell
+        (50.025, 75.0),  # Should be in a sliver cell
+    ]
+
+    found_count = 0
+    for x, y in test_points:
+        result = index.query_point(x, y, k=20)  # Use k=20 to be thorough
+        if result is not None:
+            # Verify the point is actually in the found cell
+            xv, yv, _ = grid.xyzvertices
+            verts = np.column_stack([xv[result], yv[result]])
+            from matplotlib.path import Path
+
+            path = Path(verts)
+            is_inside = path.contains_point((x, y), radius=1e-9)
+
+            if is_inside:
+                found_count += 1
+                print(f"✅ Point ({x}, {y}) correctly found in cell {result}")
+            else:
+                print(
+                    f"⚠️  Point ({x}, {y}) found in cell {result} "
+                    f"but verification failed"
+                )
+        else:
+            print(f"❌ Point ({x}, {y}) NOT FOUND")
+
+    # At least some points should be found
+    # (not all may be in cells due to Delaunay triangulation specifics)
+    assert found_count > 0, (
+        f"Should find at least some points in sliver cells, found {found_count}/3"
+    )
+
+
+def test_boundary_points():
+    """Test points exactly on cell boundaries."""
+    # Create simple 2x2 grid of triangles
+    vertices = [
+        [0, 0.0, 0.0],
+        [1, 1.0, 0.0],
+        [2, 2.0, 0.0],
+        [3, 0.0, 1.0],
+        [4, 1.0, 1.0],
+        [5, 2.0, 1.0],
+    ]
+
+    # Create triangular cells
+    cell2d = [
+        [0, 0.33, 0.33, 3, 0, 1, 3],  # Lower-left triangle
+        [1, 0.67, 0.33, 3, 1, 4, 3],  # Lower-middle triangle
+        [2, 1.33, 0.33, 3, 1, 2, 4],  # Lower-right triangle
+        [3, 1.67, 0.67, 3, 2, 5, 4],  # Upper-right triangle
+    ]
+
+    grid = VertexGrid(
+        vertices=vertices, cell2d=cell2d, top=np.ones(4) * 10.0, botm=np.zeros(4)
+    )
+
+    index = GeospatialIndex(grid)
+
+    # Test point exactly on a shared edge (between cells 0 and 1)
+    # Should find one of the two cells
+    result = index.query_point(1.0, 0.5, k=10)
+    assert result is not None or result == 0 or result == 1, (
+        "Point on boundary should be found in one of the adjacent cells"
+    )
+
+
+def test_corner_points():
+    """Test points at vertices (corners where multiple cells meet)."""
+    # Simple 2x2 grid
+    vertices = [
+        [0, 0.0, 0.0],
+        [1, 1.0, 0.0],
+        [2, 1.0, 1.0],
+        [3, 0.0, 1.0],
+        [4, 0.5, 0.5],  # Center point
+    ]
+
+    cell2d = [
+        [0, 0.5, 0.17, 3, 0, 1, 4],
+        [1, 0.83, 0.5, 3, 1, 2, 4],
+        [2, 0.5, 0.83, 3, 2, 3, 4],
+        [3, 0.17, 0.5, 3, 3, 0, 4],
+    ]
+
+    grid = VertexGrid(
+        vertices=vertices, cell2d=cell2d, top=np.ones(4) * 10.0, botm=np.zeros(4)
+    )
+
+    index = GeospatialIndex(grid)
+
+    # Test point at center vertex (shared by all 4 cells)
+    result = index.query_point(0.5, 0.5, k=10)
+    assert result is not None, (
+        "Point at center vertex should be found in one of the cells"
+    )
+    assert result in [0, 1, 2, 3], (
+        f"Result {result} should be one of the 4 center cells"
+    )
+
+
+def test_outside_grid():
+    """Test points well outside the grid."""
+    # Simple triangle
+    vertices = [
+        [0, 0.0, 0.0],
+        [1, 10.0, 0.0],
+        [2, 5.0, 10.0],
+    ]
+
+    cell2d = [
+        [0, 5.0, 3.33, 3, 0, 1, 2],
+    ]
+
+    grid = VertexGrid(
+        vertices=vertices, cell2d=cell2d, top=np.array([10.0]), botm=np.array([0.0])
+    )
+
+    index = GeospatialIndex(grid)
+
+    # Points clearly outside
+    outside_points = [
+        (-10.0, 0.0),
+        (20.0, 0.0),
+        (5.0, 20.0),
+        (15.0, 15.0),
+    ]
+
+    for x, y in outside_points:
+        result = index.query_point(x, y, k=10)
+        assert result is None, (
+            f"Point ({x}, {y}) outside grid should return None, got {result}"
+        )
+
+
+if __name__ == "__main__":
+    print("Running edge case tests...")
+
+    print("\n1. Testing thin sliver cells...")
+    test_thin_sliver_cell()
+
+    print("\n2. Testing boundary points...")
+    test_boundary_points()
+
+    print("\n3. Testing corner points...")
+    test_corner_points()
+
+    print("\n4. Testing outside grid points...")
+    test_outside_grid()
+
+    print("\n✅ All edge case tests passed!")