From a0f8fd41cd78bfdb78f050f7edb0ca2419a77bef Mon Sep 17 00:00:00 2001 From: Greg Lucas Date: Sat, 4 Jul 2026 09:04:32 -0600 Subject: [PATCH] PERF: vectorize shapely_to_path coordinate extraction Rather than pulling out each sub-geometry .xy's coordinates, we can do one get_coordinates() call and fill in the codes all at once. --- lib/cartopy/mpl/path.py | 25 +++++++++++++------------ lib/cartopy/tests/mpl/test_path.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/lib/cartopy/mpl/path.py b/lib/cartopy/mpl/path.py index 9e18a9a1f..78d4d1655 100644 --- a/lib/cartopy/mpl/path.py +++ b/lib/cartopy/mpl/path.py @@ -14,6 +14,7 @@ from matplotlib.path import Path import numpy as np +import shapely import shapely.geometry as sgeom @@ -106,20 +107,20 @@ def shapely_to_path(shape): if shape.is_empty: return Path(np.empty([0, 2])) elif isinstance(shape, sgeom.LinearRing): - return Path(np.column_stack(shape.xy), closed=True) + return Path(shapely.get_coordinates(shape), closed=True) elif isinstance(shape, (sgeom.LineString, sgeom.Point)): - return Path(np.column_stack(shape.xy)) + return Path(shapely.get_coordinates(shape)) elif isinstance(shape, sgeom.Polygon): - def poly_codes(poly): - codes = np.ones(len(poly.xy[0])) * Path.LINETO - codes[0] = Path.MOVETO - codes[-1] = Path.CLOSEPOLY - return codes - vertices = np.concatenate([np.array(shape.exterior.xy)] + - [np.array(ring.xy) for ring in - shape.interiors], 1).T - codes = np.concatenate([poly_codes(shape.exterior)] + - [poly_codes(ring) for ring in shape.interiors]) + rings = [shape.exterior, *shape.interiors] + # Shapely rings include the closing duplicate point, so len(ring.coords) + # gives the number of vertices contributed by each ring. + counts = np.fromiter((len(ring.coords) for ring in rings), dtype=np.intp, + count=len(rings)) + vertices = shapely.get_coordinates(shape) + codes = np.full(len(vertices), Path.LINETO, dtype=Path.code_type) + starts = np.cumsum(counts) - counts + codes[starts] = Path.MOVETO + codes[starts + counts - 1] = Path.CLOSEPOLY return Path(vertices, codes) elif isinstance(shape, (sgeom.MultiPolygon, sgeom.GeometryCollection, sgeom.MultiLineString, sgeom.MultiPoint)): diff --git a/lib/cartopy/tests/mpl/test_path.py b/lib/cartopy/tests/mpl/test_path.py index 269579575..19de473a6 100644 --- a/lib/cartopy/tests/mpl/test_path.py +++ b/lib/cartopy/tests/mpl/test_path.py @@ -5,6 +5,7 @@ from matplotlib.path import Path import numpy as np +from numpy.testing import assert_array_equal import pytest import shapely.geometry as sgeom @@ -102,3 +103,32 @@ def test_non_polygon_path_closing(self, path): closed_path = cpath._ensure_path_closed(path) assert isinstance(closed_path, Path) assert closed_path.vertices.size == 0 + + +class Test_shapely_to_path: + def test_polygon_with_multiple_interiors(self): + exterior = sgeom.box(0, 0, 12, 12).exterior.coords + interiors = [sgeom.box(1, 1, 2, 2, ccw=False).exterior.coords, + sgeom.box(4, 4, 5, 6, ccw=False).exterior.coords, + sgeom.box(8, 8, 9, 10, ccw=False).exterior.coords] + poly = sgeom.Polygon(exterior, interiors) + + path = cpath.shapely_to_path(poly) + + rings = [poly.exterior, *poly.interiors] + expected_vertices = np.concatenate([np.asarray(ring.coords) + for ring in rings]) + assert_array_equal(path.vertices, expected_vertices) + + # Each ring is 5 vertices long (a closed box), so the codes for + # each ring should be MOVETO, LINETO, LINETO, LINETO, CLOSEPOLY. + expected_codes = np.tile( + [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY], + len(rings)) + assert_array_equal(path.codes, expected_codes) + assert path.codes.dtype == Path.code_type + + # The path should round-trip back to an equivalent geometry. + result = cpath.path_to_shapely(path) + assert isinstance(result, sgeom.Polygon) + assert result.equals(poly)