Fix MultiPolygon patch creation for inverted inner rings (#517)

Author: timtreis

src/spatialdata_plot/pl/utils.py

@@ -52,6 +52,7 @@
 from scanpy.plotting._utils import add_colors_for_categorical_sample_annotation
 from scanpy.plotting.palettes import default_20, default_28, default_102
 from scipy.spatial import ConvexHull
+from shapely.errors import GEOSException
 from skimage.color import label2rgb
 from skimage.morphology import erosion, square
 from skimage.segmentation import find_boundaries
@@ -445,8 +446,35 @@ def _as_rgba_array(x: Any) -> np.ndarray:
     else:
         outline_c = [None] * fill_c.shape[0]
 
-    # Build DataFrame of valid geometries
-    shapes_df = pd.DataFrame(shapes, copy=True)
+    if isinstance(shapes, GeoDataFrame):
+        shapes_df: GeoDataFrame | pd.DataFrame = shapes.copy()
+    else:
+        shapes_df = pd.DataFrame(shapes, copy=True)
+
+    # Robustly normalise geometries to a canonical representation.
+    # This ensures consistent exterior/interior ring orientation so that
+    # matplotlib's fill rules handle holes correctly regardless of user input.
+    if "geometry" in shapes_df.columns:
+
+        def _normalize_geom(geom: Any) -> Any:
+            if geom is None or getattr(geom, "is_empty", False):
+                return geom
+            # shapely.normalize is available in shapely>=2; fall back to geom.normalize()
+            normalize_func = getattr(shapely, "normalize", None)
+            if callable(normalize_func):
+                try:
+                    return normalize_func(geom)
+                except (GEOSException, TypeError, ValueError):
+                    return geom
+            if hasattr(geom, "normalize"):
+                try:
+                    return geom.normalize()
+                except (GEOSException, TypeError, ValueError):
+                    return geom
+            return geom
+
+        shapes_df["geometry"] = shapes_df["geometry"].apply(_normalize_geom)
+
     shapes_df = shapes_df[shapes_df["geometry"].apply(lambda geom: not geom.is_empty)]
     shapes_df = shapes_df.reset_index(drop=True)
 
@@ -1672,52 +1700,36 @@ def _validate_polygons(shapes: GeoDataFrame) -> GeoDataFrame:
     return shapes
 
 
-def _collect_polygon_rings(
-    geom: shapely.Polygon | shapely.MultiPolygon,
-) -> list[tuple[np.ndarray, list[np.ndarray]]]:
-    """Collect exterior/interior coordinate rings from (Multi)Polygons."""
-    polygons: list[tuple[np.ndarray, list[np.ndarray]]] = []
-
-    def _collect(part: shapely.Polygon | shapely.MultiPolygon) -> None:
-        if part.geom_type == "Polygon":
-            exterior = np.asarray(part.exterior.coords)
-            interiors = [np.asarray(interior.coords) for interior in part.interiors]
-            polygons.append((exterior, interiors))
-        elif part.geom_type == "MultiPolygon":
-            for child in part.geoms:
-                _collect(child)
-        else:
-            raise ValueError(f"Unhandled geometry type: {repr(part.geom_type)}")
-
-    _collect(geom)
-    return polygons
-
+def _make_patch_from_multipolygon(mp: shapely.MultiPolygon) -> list[mpatches.PathPatch]:
+    """
+    Create PathPatches from a MultiPolygon, preserving holes robustly.
 
-def _create_ring_codes(length: int) -> npt.NDArray[np.uint8]:
-    codes = np.full(length, mpath.Path.LINETO, dtype=mpath.Path.code_type)
-    codes[0] = mpath.Path.MOVETO
-    return codes
+    This follows the same strategy as GeoPandas' internal Polygon plotting:
+    each (multi)polygon part becomes a compound Path composed of the exterior
+    ring and all interior rings. Orientation is handled by prior geometry
+    normalization rather than manual ring reversal.
+    """
+    patches: list[mpatches.PathPatch] = []
 
+    for poly in mp.geoms:
+        if poly.is_empty:
+            continue
 
-def _make_patch_from_multipolygon(mp: shapely.MultiPolygon) -> mpatches.PathPatch:
-    # https://matplotlib.org/stable/gallery/shapes_and_collections/donut.html
+        # Ensure 2D vertices in case geometries carry Z
+        exterior = np.asarray(poly.exterior.coords)[..., :2]
+        interiors = [np.asarray(ring.coords)[..., :2] for ring in poly.interiors]
 
-    patches = []
-    for exterior, interiors in _collect_polygon_rings(mp):
         if len(interiors) == 0:
+            # Simple polygon without holes
             patches.append(mpatches.Polygon(exterior, closed=True))
             continue
 
-        ring_vertices = [exterior]
-        ring_codes = [_create_ring_codes(len(exterior))]
-        for hole in interiors:
-            reversed_hole = hole[::-1]
-            ring_vertices.append(reversed_hole)
-            ring_codes.append(_create_ring_codes(len(reversed_hole)))
-
-        vertices = np.concatenate(ring_vertices)
-        all_codes = np.concatenate(ring_codes)
-        patches.append(mpatches.PathPatch(mpath.Path(vertices, all_codes)))
+        # Build a compound path: exterior + all interior rings
+        compound_path = mpath.Path.make_compound_path(
+            mpath.Path(exterior, closed=True),
+            *[mpath.Path(ring, closed=True) for ring in interiors],
+        )
+        patches.append(mpatches.PathPatch(compound_path))
 
     return patches
 

tests/_images/Shapes_can_render_multipolygon_with_inverted_inner_ring_and_disjoint_part.png

@@ -130,6 +130,81 @@ def test_plot_can_render_multipolygons_that_say_they_are_polygons(self):
 
         fig.tight_layout()
 
+    def test_plot_can_render_polygon_with_inverted_inner_ring(self):
+        ext = [
+            (7.866043666934409, 32.80184055229537),
+            (19.016191271980425, 203.48380872801957),
+            (75.90086964475744, 236.02570144190528),
+            (229.48380872801957, 235.98380872801957),
+            (235.98380872801957, 5.516191271980426),
+            (197.42585593903195, 6.144892860751103),
+            (116.5, 96.4575926540027),
+            (55.65582863082729, 12.531294107459374),
+            (7.866043666934409, 32.80184055229537),
+        ]
+
+        interior = [
+            (160.12353079731844, 173.21221665537414),
+            (181.80184055229537, 159.13395633306558),
+            (198.86604366693442, 179.80184055229537),
+            (178.19815944770465, 198.86604366693442),
+            (160.12353079731844, 173.21221665537414),
+        ]
+
+        polygon = Polygon(ext, [interior])
+        geo_df = gpd.GeoDataFrame(geometry=[polygon])
+        sdata = SpatialData(shapes={"inverted_ring": ShapesModel.parse(geo_df)})
+
+        fig, ax = plt.subplots()
+        sdata.pl.render_shapes(element="inverted_ring").pl.show(ax=ax)
+        ax.set_xlim(0, 250)
+        ax.set_ylim(0, 250)
+
+        fig.tight_layout()
+
+    def test_plot_can_render_multipolygon_with_inverted_inner_ring_and_disjoint_part(self):
+        ext = [
+            (7.866043666934409, 32.80184055229537),
+            (19.016191271980425, 203.48380872801957),
+            (75.90086964475744, 236.02570144190528),
+            (229.48380872801957, 235.98380872801957),
+            (235.98380872801957, 5.516191271980426),
+            (197.42585593903195, 6.144892860751103),
+            (116.5, 96.4575926540027),
+            (55.65582863082729, 12.531294107459374),
+            (7.866043666934409, 32.80184055229537),
+        ]
+
+        interior = [
+            (160.12353079731844, 173.21221665537414),
+            (181.80184055229537, 159.13395633306558),
+            (198.86604366693442, 179.80184055229537),
+            (178.19815944770465, 198.86604366693442),
+            (160.12353079731844, 173.21221665537414),
+        ]
+
+        # Part with a hole and non-standard orientation, plus a disjoint simple part
+        poly_with_hole = Polygon(ext, [interior])
+        disjoint_poly = Polygon(
+            [
+                (300.0, 300.0),
+                (320.0, 300.0),
+                (320.0, 320.0),
+                (300.0, 320.0),
+                (300.0, 300.0),
+            ]
+        )
+        multipoly = MultiPolygon([poly_with_hole, disjoint_poly])
+        geo_df = gpd.GeoDataFrame(geometry=[multipoly])
+        sdata = SpatialData(shapes={"inverted_ring_multipoly": ShapesModel.parse(geo_df)})
+
+        fig, ax = plt.subplots()
+        sdata.pl.render_shapes(element="inverted_ring_multipoly").pl.show(ax=ax)
+        ax.set_xlim(0, 350)
+        ax.set_ylim(0, 350)
+
+        fig.tight_layout()
+
     def test_plot_can_color_multipolygons_with_multiple_holes(self):
         square = [(0.0, 0.0), (5.0, 0.0), (5.0, 5.0), (0.0, 5.0), (0.0, 0.0)]
         first_hole = [(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)]