Skip to content

Commit

Permalink
fix(export_contours/f): support matplotlib 3.8+
Browse files Browse the repository at this point in the history 
* ContourSet.get_paths() returns one path per level
* path may contain multiple disconnected components
* walk contour components to determine connectivity
  • Loading branch information
@wpbonelli
wpbonelli committed Sep 19, 2023
1 parent 2fa32da commit c3c95ea
Show file tree
Hide file tree
Showing 2 changed files with 162 additions and 43 deletions.
15 changes: 14 additions & 1 deletion autotest/test_export.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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")
Expand Down
190 changes: 148 additions & 42 deletions flopy/export/utils.py
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,11 @@
import os
from importlib.metadata import version
from itertools import repeat
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
Expand Down Expand Up @@ -1680,6 +1683,7 @@ def export_contours(
filename: Union[str, os.PathLike],
contours,
fieldname="level",
verbose=False,
**kwargs,
):
"""
Expand All @@ -1693,6 +1697,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
Expand All @@ -1706,20 +1712,66 @@ def export_contours(
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(
Expand Down Expand Up @@ -1759,54 +1811,108 @@ def export_contourf(
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(
Expand Down

0 comments on commit c3c95ea

Please sign in to comment.