update export_contours (can't rely on a path for each disjoint component)

Author: wpbonelli

autotest/test_export.py

@@ -51,7 +51,6 @@
 from flopy.utils.crs import get_authority_crs
 from flopy.utils.geometry import Polygon
 
-
 HAS_PYPROJ = has_pkg("pyproj", strict=True)
 if HAS_PYPROJ:
     import pyproj
@@ -717,14 +716,14 @@ def test_export_contours(function_tmpdir, example_data_path):
     with Reader(filename) as r:
         shapes = r.shapes()
         # expect 65 with standard mpl contours (structured grids), 86 with tricontours
-        # assert len(shapes) >= 65
+        assert len(shapes) >= 65
 
         # debugging
-        for s in shapes:
-            x = [i[0] for i in s.points[:]]
-            y = [i[1] for i in s.points[:]]
-            plt.plot(x, y)
-        plt.show()
+        # for s in shapes:
+        #     x = [i[0] for i in s.points[:]]
+        #     y = [i[1] for i in s.points[:]]
+        #     plt.plot(x, y)
+        # plt.show()
 
 
 @pytest.mark.mf6

flopy/export/utils.py

@@ -1,8 +1,9 @@
 import os
 from pathlib import Path
-from typing import Optional, Union
+from typing import Union
 
 import numpy as np
+from matplotlib.path import Path
 
 from ..datbase import DataInterface, DataListInterface, DataType
 from ..mbase import BaseModel, ModelInterface
@@ -1706,20 +1707,55 @@ def export_contours(
     if not isinstance(contours, list):
         contours = [contours]
 
+    # Below we export a geometry (shape) for each
+    # disjoint contour component. Multiple shapes
+    # may correspond to the same contour level.
     geoms = []
     level = []
+
     for ctr in contours:
-        levels = ctr.levels
-        for i, c in enumerate(ctr.collections):
-            paths = c.get_paths()
-            geoms += [LineString(p.vertices) for p in paths]
-            level += list(np.ones(len(paths)) * levels[i])
+        # ContourSet.collections was deprecated with
+        # matplotlib 3.8, ContourSet is a collection
+        # of Paths now. Each Path corresponds to one
+        # contour level, and may contain one or more
+        # (possibly disconnected) components.
+
+        paths = ctr.get_paths()
+        for pi, path in enumerate(paths):
+            # skip empty paths
+            if path.vertices.shape[0] == 0:
+                continue
+
+            # Each Path may contain multiple components
+            # so we unpack them as separate geometries.
+            comps = []
+            segs = []
+            for seg in path.iter_segments():
+                pts, code = seg
+                if code == Path.MOVETO:
+                    if len(segs) > 0:
+                        comps.append(LineString(segs))
+                        segs = []
+                    segs.append(pts)
+                elif code == Path.LINETO:
+                    segs.append(pts)
+                elif code == Path.CLOSEPOLY:
+                    segs.append(pts)
+                    segs.append(segs[0])  # add closing segment
+                    comps.append(LineString(segs))
+                    segs = []
+            if len(segs) > 0:
+                comps.append(LineString(segs))
+            if len(segs) == 1:
+                raise RuntimeError(
+                    f"Encountered component with only 1 vertex. This should not have happened, contact developers."
+                )
+            level += list(np.ones(len(comps)) * ctr.levels[pi])
+            geoms += comps
 
-    # convert the dictionary to a recarray
     ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
 
     recarray2shp(ra, geoms, filename, **kwargs)
-    return
 
 
 def export_contourf(
@@ -1806,7 +1842,6 @@ def export_contourf(
     ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
 
     recarray2shp(ra, geoms, filename, **kwargs)
-    return
 
 
 def export_array_contours(