From cd4f58d3233c206ed98f1591e8bc6d740b3bcc42 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Sat, 16 Sep 2023 14:36:57 -0400 Subject: [PATCH 1/3] fix(export_contours/f): support matplotlib 3.8+ * ContourSet.get_paths() returns one path per level * path may contain multiple disconnected components * walk contour components to determine connectivity --- autotest/test_export.py | 15 ++- flopy/export/utils.py | 196 +++++++++++++++++++++++++++++++--------- 2 files changed, 168 insertions(+), 43 deletions(-) diff --git a/autotest/test_export.py b/autotest/test_export.py index afe9bdf39a..c7f84661da 100644 --- a/autotest/test_export.py +++ b/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") diff --git a/flopy/export/utils.py b/flopy/export/utils.py index e8266787a4..2d90a4f30e 100644 --- a/flopy/export/utils.py +++ b/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,6 +1695,8 @@ 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 @@ -1700,26 +1704,76 @@ def export_contours( 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,115 @@ 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 + + poly = None + polys = [] + + def append_poly(segs): + # TODO: handle holes? + # check if this is a multipolygon by checking vertex + # order. + # if is_clockwise(segs): + # # Clockwise is a hole, set to interiors + # if not poly.interiors: + # poly.interiors = [segs] + # else: + # poly.interiors.append(segs) + # else: + poly = Polygon(segs) + polys.append(poly) + + segs = [] + for seg in path.iter_segments(): + pts, code = seg + if code == Path.MOVETO: + if len(segs) > 0: + append_poly(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 + append_poly(segs) + segs = [] + if len(segs) > 0: + append_poly(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( From 0a5b89e7c6bc28aeb0ea8c31cc1271907faf5ab3 Mon Sep 17 00:00:00 2001 From: Wes Bonelli Date: Tue, 19 Sep 2023 09:37:24 -0400 Subject: [PATCH 2/3] remove commented interior hole handling for mpl 3.8+ --- flopy/export/utils.py | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 2d90a4f30e..879c4fde5a 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -1875,29 +1875,13 @@ def export_contourf( if path.vertices.shape[0] == 0: continue - poly = None polys = [] - - def append_poly(segs): - # TODO: handle holes? - # check if this is a multipolygon by checking vertex - # order. - # if is_clockwise(segs): - # # Clockwise is a hole, set to interiors - # if not poly.interiors: - # poly.interiors = [segs] - # else: - # poly.interiors.append(segs) - # else: - poly = Polygon(segs) - polys.append(poly) - segs = [] for seg in path.iter_segments(): pts, code = seg if code == Path.MOVETO: if len(segs) > 0: - append_poly(segs) + polys.append(Polygon(segs)) segs = [] segs.append(pts) elif code == Path.LINETO: @@ -1905,10 +1889,10 @@ def append_poly(segs): elif code == Path.CLOSEPOLY: segs.append(pts) segs.append(segs[0]) # add closing segment - append_poly(segs) + polys.append(Polygon(segs)) segs = [] if len(segs) > 0: - append_poly(segs) + polys.append(Polygon(segs)) geoms.extend(polys) level.extend(repeat(ctr.levels[pi], len(polys))) From 5a41f5cdd11b3dedb0f2286d378dad3d2fda1ffc Mon Sep 17 00:00:00 2001 From: Wes Bonelli Date: Tue, 19 Sep 2023 10:16:21 -0400 Subject: [PATCH 3/3] fix mpl version comparisons --- flopy/export/utils.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 879c4fde5a..0277e55d43 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -4,6 +4,7 @@ from typing import Union import numpy as np +from packaging.version import Version from ..datbase import DataInterface, DataListInterface, DataType from ..mbase import BaseModel, ModelInterface @@ -1726,12 +1727,12 @@ def export_contours( # (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 + # but post-3.8, we have to walk the segments # to distinguish disconnected components. - mpl_ver = version("matplotlib") + mpl_ver = Version(version("matplotlib")) for ctr in contours: - if mpl_ver < "3.8.0": + if mpl_ver < Version("3.8.0"): levels = ctr.levels for i, c in enumerate(ctr.collections): paths = c.get_paths() @@ -1831,12 +1832,12 @@ def export_contourf( # (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 + # but post-3.8, we have to walk the segments # to distinguish disconnected components. - mpl_ver = version("matplotlib") + mpl_ver = Version(version("matplotlib")) for ctr in contours: - if mpl_ver < "3.8.0": + if mpl_ver < Version("3.8.0"): levels = ctr.levels for idx, col in enumerate(ctr.collections): for contour_path in col.get_paths():