Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
92 changes: 52 additions & 40 deletions src/spatialdata_plot/pl/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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

Expand Down
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
75 changes: 75 additions & 0 deletions tests/pl/test_render_shapes.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)]
Expand Down