fix(export_contours/f): support matplotlib 3.8+

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
@@ -683,6 +682,13 @@ def test_export_contourf(function_tmpdir, example_data_path):
             len(shapes) >= 65
         ), "multipolygons were skipped in contourf routine"
 
+        # 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()
+
 
 @pytest.mark.mf6
 @requires_pkg("shapefile", "shapely")
@@ -712,6 +718,13 @@ def test_export_contours(function_tmpdir, example_data_path):
         # expect 65 with standard mpl contours (structured grids), 86 with tricontours
         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()
+
 
 @pytest.mark.mf6
 @requires_pkg("shapely")

flopy/export/utils.py

@@ -1,6 +1,7 @@
 import os
+from itertools import repeat
 from pathlib import Path
-from typing import Optional, Union
+from typing import Union
 
 import numpy as np
 
@@ -1680,6 +1681,7 @@ def export_contours(
     filename: Union[str, os.PathLike],
     contours,
     fieldname="level",
+    verbose=False,
     **kwargs,
 ):
     """
@@ -1693,33 +1695,85 @@ def export_contours(
         (object returned by matplotlib.pyplot.contour)
     fieldname : str
         gis attribute table field name
+    verbose : bool, optional, default False
+        whether to show verbose output
     **kwargs : key-word arguments to flopy.export.shapefile_utils.recarray2shp
 
     Returns
     -------
     df : dataframe of shapefile contents
 
     """
+    from importlib.metadata import version
+
+    from matplotlib.path import Path
+
     from ..utils.geometry import LineString
     from .shapefile_utils import recarray2shp
 
     if not isinstance(contours, list):
         contours = [contours]
 
+    # Export a linestring for each contour component.
+    # Levels may have multiple disconnected components.
     geoms = []
     level = []
+
+    # ContourSet.collections was deprecated with
+    # matplotlib 3.8. ContourSet is a collection
+    # of Paths, where each Path corresponds to a
+    # contour level, and may contain one or more
+    # (possibly disconnected) components. Before
+    # 3.8, iterating over ContourSet.collections
+    # and enumerating from get_paths() suffices,
+    # but post-3.8 we have to walk the segmments
+    # to distinguish disconnected components.
+    mpl_ver = version("matplotlib")
+
     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])
+        if mpl_ver < "3.8.0":
+            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])
+        else:
+            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.
+                lines = []
+                segs = []
+                for seg in path.iter_segments():
+                    pts, code = seg
+                    if code == Path.MOVETO:
+                        if len(segs) > 0:
+                            lines.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
+                        lines.append(LineString(segs))
+                        segs = []
+                if len(segs) > 0:
+                    lines.append(LineString(segs))
+
+                level += list(np.ones(len(lines)) * ctr.levels[pi])
+                geoms += lines
+
+    if verbose:
+        print(f"Writing {len(level)} contour lines")
 
-    # 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(
@@ -1756,57 +1810,99 @@ def export_contourf(
     >>> export_contourf('myfilledcontours.shp', cs)
 
     """
+    from importlib.metadata import version
+
+    from matplotlib.path import Path
+
     from ..utils.geometry import Polygon, is_clockwise
     from .shapefile_utils import recarray2shp
 
+    if not isinstance(contours, list):
+        contours = [contours]
+
+    # export a polygon for each filled contour
     geoms = []
     level = []
 
-    if not isinstance(contours, list):
-        contours = [contours]
+    # ContourSet.collections was deprecated with
+    # matplotlib 3.8. ContourSet is a collection
+    # of Paths, where each Path corresponds to a
+    # contour level, and may contain one or more
+    # (possibly disconnected) components. Before
+    # 3.8, iterating over ContourSet.collections
+    # and enumerating from get_paths() suffices,
+    # but post-3.8 we have to walk the segmments
+    # to distinguish disconnected components.
+    mpl_ver = version("matplotlib")
 
-    for c in contours:
-        levels = c.levels
-        for idx, col in enumerate(c.collections):
-            # Loop through all polygons that have the same intensity level
-            for contour_path in col.get_paths():
-                # Create the polygon(s) for this intensity level
-                poly = None
-                for ncp, cp in enumerate(contour_path.to_polygons()):
-                    x = cp[:, 0]
-                    y = cp[:, 1]
-                    verts = [(i[0], i[1]) for i in zip(x, y)]
-                    new_shape = Polygon(verts)
-
-                    if ncp == 0:
-                        poly = new_shape
-                    else:
-                        # check if this is a multipolygon by checking vertex
-                        # order.
-                        if is_clockwise(verts):
-                            # Clockwise is a hole, set to interiors
-                            if not poly.interiors:
-                                poly.interiors = [new_shape.exterior]
-                            else:
-                                poly.interiors.append(new_shape.exterior)
-                        else:
-                            geoms.append(poly)
-                            level.append(levels[idx])
+    for ctr in contours:
+        if mpl_ver < "3.8.0":
+            levels = ctr.levels
+            for idx, col in enumerate(ctr.collections):
+                for contour_path in col.get_paths():
+                    # Create the polygon(s) for this intensity level
+                    poly = None
+                    for ncp, cp in enumerate(contour_path.to_polygons()):
+                        x = cp[:, 0]
+                        y = cp[:, 1]
+                        verts = [(i[0], i[1]) for i in zip(x, y)]
+                        new_shape = Polygon(verts)
+
+                        if ncp == 0:
                             poly = new_shape
+                        else:
+                            # check if this is a multipolygon by checking vertex
+                            # order.
+                            if is_clockwise(verts):
+                                # Clockwise is a hole, set to interiors
+                                if not poly.interiors:
+                                    poly.interiors = [new_shape.exterior]
+                                else:
+                                    poly.interiors.append(new_shape.exterior)
+                            else:
+                                geoms.append(poly)
+                                level.append(levels[idx])
+                                poly = new_shape
 
-                if poly is not None:
-                    # store geometry object
-                    geoms.append(poly)
+                    if poly is not None:
+                        # store geometry object
+                        geoms.append(poly)
                     level.append(levels[idx])
+        else:
+            paths = ctr.get_paths()
+            for pi, path in enumerate(paths):
+                # skip empty paths
+                if path.vertices.shape[0] == 0:
+                    continue
+
+                polys = []
+                segs = []
+                for seg in path.iter_segments():
+                    pts, code = seg
+                    if code == Path.MOVETO:
+                        if len(segs) > 0:
+                            polys.append(Polygon(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
+                        polys.append(Polygon(segs))
+                        segs = []
+                if len(segs) > 0:
+                    polys.append(Polygon(segs))
+
+                geoms.extend(polys)
+                level.extend(repeat(ctr.levels[pi], len(polys)))
 
     if verbose:
         print(f"Writing {len(level)} polygons")
 
-    # Create recarray
     ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray)
 
     recarray2shp(ra, geoms, filename, **kwargs)
-    return
 
 
 def export_array_contours(