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 cd4f58d
Show file tree
Hide file tree
Showing 2 changed files with 168 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
196 changes: 154 additions & 42 deletions flopy/export/utils.py
View file
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -1680,6 +1681,7 @@ def export_contours(
filename: Union[str, os.PathLike],
contours,
fieldname="level",
verbose=False,
**kwargs,
):
"""
Expand All @@ -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(
Expand Down Expand Up @@ -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(
Expand Down

0 comments on commit cd4f58d

Please sign in to comment.