From 5f60c0997a8862fc69f16ac46cfc24428e62e116 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Mon, 15 Apr 2024 19:35:56 +0200
Subject: [PATCH] fix(#2152): improve gridintersect geometry creation for
 vertex grids - improve performance _vtx_grid_to_geoms_cellids()

---
 flopy/utils/gridintersect.py | 71 +++++++++++-------------------------
 1 file changed, 21 insertions(+), 50 deletions(-)

diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index 8b3c5cda2f..064f162aec 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -429,57 +429,28 @@ def _vtx_grid_to_geoms_cellids(self):
         cellids : array_like
             array of cellids
         """
-        shapely_geo = import_optional_dependency("shapely.geometry")
-
-        # for cell2d rec-arrays
-        geoms = []
-        cellids = []
-        if isinstance(self.mfgrid._cell2d, np.recarray):
-            for icell in self.mfgrid._cell2d.icell2d:
-                points = []
-                icverts = [
-                    f"icvert_{i}"
-                    for i in range(self.mfgrid._cell2d["ncvert"][icell])
-                ]
-                for iv in self.mfgrid._cell2d[icverts][icell]:
-                    if self.local:
-                        xy = (
-                            self.mfgrid._vertices.xv[iv],
-                            self.mfgrid._vertices.yv[iv],
-                        )
-                    else:
-                        xy = (
-                            self.mfgrid.verts[iv, 0],
-                            self.mfgrid.verts[iv, 1],
-                        )
-                    points.append(xy)
-                # close the polygon, if necessary
-                if points[0] != points[-1]:
-                    points.append(points[0])
-                geoms.append(shapely_geo.Polygon(points))
-                cellids.append(icell)
-        # for cell2d lists
-        elif isinstance(self.mfgrid._cell2d, list):
-            for icell in range(len(self.mfgrid._cell2d)):
-                points = []
-                for iv in self.mfgrid._cell2d[icell][4:]:
-                    if self.local:
-                        xy = (
-                            self.mfgrid._vertices[iv][1],
-                            self.mfgrid._vertices[iv][2],
-                        )
-                    else:
-                        xy = (
-                            self.mfgrid.verts[iv, 0],
-                            self.mfgrid.verts[iv, 1],
+        shapely = import_optional_dependency("shapely")
+        if self.local:
+            geoms = [
+                shapely.polygons(
+                    list(
+                        zip(
+                            *self.mfgrid.get_local_coords(
+                                *np.array(
+                                    self.mfgrid.get_cell_vertices(node)
+                                ).T
+                            )
                         )
-                    points.append(xy)
-                # close the polygon, if necessary
-                if points[0] != points[-1]:
-                    points.append(points[0])
-                geoms.append(shapely_geo.Polygon(points))
-                cellids.append(icell)
-        return np.array(geoms), np.array(cellids)
+                    )
+                )
+                for node in range(self.mfgrid.nnodes)
+            ]
+        else:
+            geoms = [
+                shapely.polygons(self.mfgrid.get_cell_vertices(node))
+                for node in range(self.mfgrid.nnodes)
+            ]
+        return np.array(geoms), np.arange(self.mfgrid.nnodes)
 
     def _rect_grid_to_shape_list(self):
         """internal method, list of shapely polygons for structured grid cells.