diff --git a/src/spatialdata_plot/pl/utils.py b/src/spatialdata_plot/pl/utils.py index cf21d212..cccb587e 100644 --- a/src/spatialdata_plot/pl/utils.py +++ b/src/spatialdata_plot/pl/utils.py @@ -1246,47 +1246,52 @@ def _get_linear_colormap(colors: list[str], background: str) -> list[LinearSegme return [LinearSegmentedColormap.from_list(c, [background, c], N=256) for c in colors] -def _split_multipolygon_into_outer_and_inner(mp: shapely.MultiPolygon): # type: ignore - # https://stackoverflow.com/a/21922058 - - for geom in mp.geoms: - if geom.geom_type == "MultiPolygon": - exterior_coords = [] - interior_coords = [] - for part in geom: - epc = _split_multipolygon_into_outer_and_inner(part) # Recursive call - exterior_coords += epc["exterior_coords"] - interior_coords += epc["interior_coords"] - elif geom.geom_type == "Polygon": - exterior_coords = geom.exterior.coords[:] - interior_coords = [] - for interior in geom.interiors: - interior_coords += interior.coords[:] +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(geom.type)}") + raise ValueError(f"Unhandled geometry type: {repr(part.geom_type)}") - return interior_coords, exterior_coords + _collect(geom) + return polygons + + +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 def _make_patch_from_multipolygon(mp: shapely.MultiPolygon) -> mpatches.PathPatch: # https://matplotlib.org/stable/gallery/shapes_and_collections/donut.html patches = [] - for geom in mp.geoms: - if len(geom.interiors) == 0: - # polygon has no holes - patches += [mpatches.Polygon(geom.exterior.coords, closed=True)] - else: - inside, outside = _split_multipolygon_into_outer_and_inner(mp) - if len(inside) > 0: - codes = np.ones(len(inside), dtype=mpath.Path.code_type) * mpath.Path.LINETO - codes[0] = mpath.Path.MOVETO - all_codes = np.concatenate((codes, codes)) - vertices = np.concatenate((outside, inside[::-1])) - else: - all_codes = [] - vertices = np.concatenate(outside) - patches += [mpatches.PathPatch(mpath.Path(vertices, all_codes))] + for exterior, interiors in _collect_polygon_rings(mp): + if len(interiors) == 0: + 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))) return patches diff --git a/tests/_images/Shapes_can_color_multipolygons_with_multiple_holes.png b/tests/_images/Shapes_can_color_multipolygons_with_multiple_holes.png new file mode 100644 index 00000000..39d3163a Binary files /dev/null and b/tests/_images/Shapes_can_color_multipolygons_with_multiple_holes.png differ diff --git a/tests/_images/Shapes_can_render_multipolygons_with_multiple_holes.png b/tests/_images/Shapes_can_render_multipolygons_with_multiple_holes.png new file mode 100644 index 00000000..14f75db3 Binary files /dev/null and b/tests/_images/Shapes_can_render_multipolygons_with_multiple_holes.png differ diff --git a/tests/pl/test_render_shapes.py b/tests/pl/test_render_shapes.py index e8b301a1..11984397 100644 --- a/tests/pl/test_render_shapes.py +++ b/tests/pl/test_render_shapes.py @@ -99,6 +99,46 @@ def _make_multi(): sdata["table"] = table sdata.pl.render_shapes(color="val", fill_alpha=0.3).pl.show() + def test_plot_can_render_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)] + second_hole = [(3.0, 3.0), (4.0, 3.0), (4.0, 4.0), (3.0, 4.0), (3.0, 3.0)] + multipoly = MultiPolygon([Polygon(square, holes=[first_hole, second_hole])]) + cell_polygon_table = gpd.GeoDataFrame(geometry=gpd.GeoSeries([multipoly])) + sd_polygons = ShapesModel.parse(cell_polygon_table) + sdata = SpatialData(shapes={"two_holes": sd_polygons}) + + fig, ax = plt.subplots() + sdata.pl.render_shapes(element="two_holes").pl.show(ax=ax) + ax.set_xlim(-1, 6) + ax.set_ylim(-1, 6) + + 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)] + second_hole = [(3.0, 3.0), (4.0, 3.0), (4.0, 4.0), (3.0, 4.0), (3.0, 3.0)] + multipoly = MultiPolygon([Polygon(square, holes=[first_hole, second_hole])]) + cell_polygon_table = gpd.GeoDataFrame(geometry=gpd.GeoSeries([multipoly])) + cell_polygon_table["instance_id"] = [0] + sd_polygons = ShapesModel.parse(cell_polygon_table) + + adata = anndata.AnnData(pd.DataFrame({"value": [1]})) + adata.obs["region"] = pd.Categorical(["two_holes"] * adata.n_obs) + adata.obs["instance_id"] = [0] + adata.obs["category"] = ["holey"] + table = TableModel.parse(adata, region="two_holes", region_key="region", instance_key="instance_id") + + sdata = SpatialData(shapes={"two_holes": sd_polygons}, tables={"table": table}) + + fig, ax = plt.subplots() + sdata.pl.render_shapes(element="two_holes", color="category", table_name="table").pl.show(ax=ax) + ax.set_xlim(-1, 6) + ax.set_ylim(-1, 6) + + fig.tight_layout() + def test_plot_can_color_from_geodataframe(self, sdata_blobs: SpatialData): blob = deepcopy(sdata_blobs) blob["table"].obs["region"] = pd.Categorical(["blobs_polygons"] * blob["table"].n_obs)