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.

Author: greglucas

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)):

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)