remove structured

Author: adacovsk

flopy/utils/geospatial_index.py

@@ -1,9 +1,12 @@
 """
-Unified geospatial indexing for FloPy grids.
+Geospatial indexing for FloPy vertex and unstructured grids.
 
-Provides efficient spatial queries for all grid types using KD-tree
-with cell centroids AND vertices for robust edge case handling,
-plus pre-computed ConvexHull equations for fast point-in-polygon testing.
+Provides efficient spatial queries using KD-tree with cell centroids
+AND vertices for robust edge case handling, plus pre-computed ConvexHull
+equations for fast point-in-polygon testing.
+
+Note: StructuredGrid has its own optimized spatial methods and does not
+use this index.
 """
 
 import numpy as np
@@ -12,7 +15,7 @@
 
 class GeospatialIndex:
     """
-    Geospatial index for efficient geometric queries on model grids.
+    Geospatial index for efficient geometric queries on vertex/unstructured grids.
 
     Uses KD-tree indexing with cell centroids + vertices to find candidate cells,
     then pre-computed ConvexHull hyperplane equations for fast vectorized
@@ -22,10 +25,13 @@ class GeospatialIndex:
     - Points near cell boundaries
     - Points in thin/sliver cells where centroid is far from the cell
 
+    Note: StructuredGrid has its own optimized spatial methods and should not
+    use this index.
+
     Parameters
     ----------
-    grid : Grid
-        Any FloPy grid object (StructuredGrid, VertexGrid, UnstructuredGrid)
+    grid : VertexGrid or UnstructuredGrid
+        A FloPy vertex or unstructured grid object
     epsilon : float, optional
         Tolerance for point-in-cell tests. Used for both bounding box
         expansion and ConvexHull hyperplane distance tests. Ensures boundary
@@ -52,10 +58,11 @@ class GeospatialIndex:
 
     Examples
     --------
-    >>> from flopy.discretization import StructuredGrid
+    >>> from flopy.discretization import VertexGrid
     >>> from flopy.utils.geospatial_index import GeospatialIndex
     >>>
-    >>> grid = StructuredGrid(delr=np.ones(10), delc=np.ones(10))
+    >>> # Create a simple triangular grid
+    >>> grid = VertexGrid(vertices, cell2d)
     >>> index = GeospatialIndex(grid)
     >>>
     >>> # Single point query
@@ -67,12 +74,12 @@ class GeospatialIndex:
 
     def __init__(self, grid, epsilon=1e-3):
         """
-        Build geospatial index for a model grid.
+        Build geospatial index for a vertex or unstructured grid.
 
         Parameters
         ----------
-        grid : Grid
-            Any FloPy grid (StructuredGrid, VertexGrid, UnstructuredGrid)
+        grid : VertexGrid or UnstructuredGrid
+            A FloPy vertex or unstructured grid
         epsilon : float, optional
             Tolerance for point-in-cell tests.
             Used for bounding box expansion and ConvexHull tests.
@@ -97,37 +104,37 @@ def _build_index(self):
         points = []
         point_to_cell = []
 
-        # Get grid dimensions
-        if self.grid.grid_type == "structured":
-            self.ncells = self.grid.nrow * self.grid.ncol
-            self._build_structured_index(points, point_to_cell)
-        elif self.grid.grid_type in ("vertex", "unstructured"):
-            if self.is_3d:
-                # 3D indexing: index all nnodes with x,y,z coordinates
-                self.ncells = self.grid.nnodes
-                self._build_3d_index(points, point_to_cell)
-            else:
-                # 2D indexing: index only 2D cells with x,y coordinates
-                if hasattr(self.grid, "ncpl"):
-                    ncpl = self.grid.ncpl
-                    if ncpl is None:
-                        # ncpl not set, fall back to xcellcenters
-                        self.ncells = len(self.grid.xcellcenters)
-                    elif np.isscalar(ncpl):
-                        # VertexGrid: ncpl is scalar
-                        self.ncells = ncpl
-                    else:
-                        # UnstructuredGrid: ncpl is array
-                        # Use first layer's cell count for 2D spatial indexing
-                        if len(ncpl) > 0:
-                            self.ncells = ncpl[0]
-                        else:
-                            self.ncells = len(self.grid.xcellcenters)
-                else:
-                    self.ncells = len(self.grid.xcellcenters)
-                self._build_vertex_index(points, point_to_cell)
+        # Get grid dimensions for vertex/unstructured grids
+        if self.grid.grid_type not in ("vertex", "unstructured"):
+            raise ValueError(
+                f"GeospatialIndex only supports vertex and unstructured grids, "
+                f"got: {self.grid.grid_type}"
+            )
+
+        if self.is_3d:
+            # 3D indexing: index all nnodes with x,y,z coordinates
+            self.ncells = self.grid.nnodes
+            self._build_3d_index(points, point_to_cell)
         else:
-            raise ValueError(f"Unknown grid type: {self.grid.grid_type}")
+            # 2D indexing: index only 2D cells with x,y coordinates
+            if hasattr(self.grid, "ncpl"):
+                ncpl = self.grid.ncpl
+                if ncpl is None:
+                    # ncpl not set, fall back to xcellcenters
+                    self.ncells = len(self.grid.xcellcenters)
+                elif np.isscalar(ncpl):
+                    # VertexGrid: ncpl is scalar
+                    self.ncells = ncpl
+                else:
+                    # UnstructuredGrid: ncpl is array
+                    # Use first layer's cell count for 2D spatial indexing
+                    if len(ncpl) > 0:
+                        self.ncells = ncpl[0]
+                    else:
+                        self.ncells = len(self.grid.xcellcenters)
+            else:
+                self.ncells = len(self.grid.xcellcenters)
+            self._build_vertex_index(points, point_to_cell)
 
         # Build KD-tree from centroids + vertices
         self.points = np.array(points)
@@ -141,26 +148,6 @@ def _build_index(self):
             # For 3D, store z-bounds for each cell
             self._precompute_3d_bounds()
 
-    def _build_structured_index(self, points, point_to_cell):
-        """Build index for structured grid - centroids + vertices."""
-        xc = self.grid.xcellcenters
-        yc = self.grid.ycellcenters
-        xv = self.grid.xvertices
-        yv = self.grid.yvertices
-
-        for i in range(self.grid.nrow):
-            for j in range(self.grid.ncol):
-                cellid = i * self.grid.ncol + j
-
-                # Add centroid
-                points.append([xc[i, j], yc[i, j]])
-                point_to_cell.append(cellid)
-
-                # Add 4 corner vertices
-                for vi, vj in [(i, j), (i, j + 1), (i + 1, j + 1), (i + 1, j)]:
-                    points.append([xv[vi, vj], yv[vi, vj]])
-                    point_to_cell.append(cellid)
-
     def _build_vertex_index(self, points, point_to_cell):
         """Build index for vertex/unstructured grid - centroids + vertices."""
         # Disable copy cache to avoid expensive deepcopy during index build
@@ -284,24 +271,8 @@ def _precompute_3d_bounds(self):
 
     def _get_cell_vertices(self, cellid):
         """Get vertices for a cell."""
-        if self.grid.grid_type == "structured":
-            i = cellid // self.grid.ncol
-            j = cellid % self.grid.ncol
-            xv = self.grid.xvertices
-            yv = self.grid.yvertices
-            return np.array(
-                [
-                    [xv[i, j], yv[i, j]],
-                    [xv[i, j + 1], yv[i, j + 1]],
-                    [xv[i + 1, j + 1], yv[i + 1, j + 1]],
-                    [xv[i + 1, j], yv[i + 1, j]],
-                ]
-            )
-        elif self.grid.grid_type in ("vertex", "unstructured"):
-            xv, yv, _ = self.grid.xyzvertices
-            return np.column_stack([xv[cellid], yv[cellid]])
-        else:
-            raise ValueError(f"Unknown grid type: {self.grid.grid_type}")
+        xv, yv, _ = self.grid.xyzvertices
+        return np.column_stack([xv[cellid], yv[cellid]])
 
     def query_point(self, x, y, z=None, k=None):
         """
@@ -540,11 +511,9 @@ def _point_in_cell_3d(self, point, cellid):
 
     def _apply_tiebreaker(self, matching_cells):
         """
-        Apply grid-specific tie-breaking when multiple cells match.
+        Apply 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
+        For vertex/unstructured grids: choose cell with lowest cell ID.
 
         Parameters
         ----------
@@ -559,21 +528,7 @@ def _apply_tiebreaker(self, matching_cells):
         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)
+        return min(matching_cells)
 
     def _find_layer_for_z(self, icell2d, z):
         """