diff --git a/src/spatialdata_plot/pl/utils.py b/src/spatialdata_plot/pl/utils.py index 954afe40..97ffac1d 100644 --- a/src/spatialdata_plot/pl/utils.py +++ b/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 diff --git a/tests/_images/Shapes_can_render_multipolygon_with_inverted_inner_ring_and_disjoint_part.png b/tests/_images/Shapes_can_render_multipolygon_with_inverted_inner_ring_and_disjoint_part.png new file mode 100644 index 00000000..1b8b3fde Binary files /dev/null and b/tests/_images/Shapes_can_render_multipolygon_with_inverted_inner_ring_and_disjoint_part.png differ diff --git a/tests/_images/Shapes_can_render_polygon_with_inverted_inner_ring.png b/tests/_images/Shapes_can_render_polygon_with_inverted_inner_ring.png new file mode 100644 index 00000000..d5bd7964 Binary files /dev/null and b/tests/_images/Shapes_can_render_polygon_with_inverted_inner_ring.png differ diff --git a/tests/pl/test_render_shapes.py b/tests/pl/test_render_shapes.py index 1acf3d3c..fce236d6 100644 --- a/tests/pl/test_render_shapes.py +++ b/tests/pl/test_render_shapes.py @@ -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)]