Support MultiPolygons in render_shapes()

CHANGELOG.md

@@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning][].
 
 ### Added
 
+-   Multipolygons are now handled correctly (#93)
+
 ### Fixed
 
 -   Legend order is now deterministic (#143)

src/spatialdata_plot/pl/render.py

@@ -6,6 +6,7 @@
 from functools import partial
 from typing import Any, Callable, Union
 
+import geopandas as gpd
 import matplotlib
 import numpy as np
 import pandas as pd
@@ -35,6 +36,7 @@
     _decorate_axs,
     _get_colors_for_categorical_obs,
     _get_linear_colormap,
+    _make_patch_from_multipolygon,
     _map_color_seg,
     _maybe_set_colors,
     _normalize,
@@ -119,18 +121,22 @@ def _get_collection_shape(
         outline_alpha: None | float = None,
         **kwargs: Any,
     ) -> PatchCollection:
-        patches = []
-        for shape in shapes:
-            # remove empty points/polygons
-            shape = shape[shape["geometry"].apply(lambda geom: not geom.is_empty)]
-            # We assume that all elements in one collection are of the same type
-            if shape["geometry"].iloc[0].geom_type == "Polygon":
-                patches += [Polygon(p.exterior.coords, closed=True) for p in shape["geometry"]]
-            elif shape["geometry"].iloc[0].geom_type == "Point":
-                patches += [
-                    Circle((circ.x, circ.y), radius=r * s) for circ, r in zip(shape["geometry"], shape["radius"])
-                ]
-
+        """
+        Get a PatchCollection for rendering given geometries with specified colors and outlines.
+
+        Args:
+        - shapes (list[GeoDataFrame]): List of geometrical shapes.
+        - c: Color parameter.
+        - s (float): Size of the shape.
+        - norm: Normalization for the color map.
+        - fill_alpha (float, optional): Opacity for the fill color.
+        - outline_alpha (float, optional): Opacity for the outline.
+        - **kwargs: Additional keyword arguments.
+
+        Returns
+        -------
+        - PatchCollection: Collection of patches for rendering.
+        """
         cmap = kwargs["cmap"]
 
         try:
@@ -149,16 +155,60 @@ def _get_collection_shape(
         if render_params.outline_params.outline:
             outline_c = ColorConverter().to_rgba_array(render_params.outline_params.outline_color)
             outline_c[..., -1] = render_params.outline_alpha
+            outline_c = outline_c.tolist()
         else:
-            outline_c = None
+            outline_c = [None]
+        outline_c = outline_c * fill_c.shape[0]
+
+        shapes_df = pd.DataFrame(shapes, copy=True)
+
+        # remove empty points/polygons
+        shapes_df = shapes_df[shapes_df["geometry"].apply(lambda geom: not geom.is_empty)]
+
+        rows = []
+
+        def assign_fill_and_outline_to_row(
+            shapes: list[GeoDataFrame], fill_c: list[Any], outline_c: list[Any], row: pd.Series, idx: int
+        ) -> None:
+            if len(shapes) > 1 and len(fill_c) == 1:
+                row["fill_c"] = fill_c
+                row["outline_c"] = outline_c
+            else:
+                row["fill_c"] = fill_c[idx]
+                row["outline_c"] = outline_c[idx]
+
+        # Match colors to the geometry, potentially expanding the row in case of
+        # multipolygons
+        for idx, row in shapes_df.iterrows():
+            geom = row["geometry"]
+            if geom.geom_type == "Polygon":
+                row = row.to_dict()
+                row["geometry"] = Polygon(geom.exterior.coords, closed=True)
+                assign_fill_and_outline_to_row(shapes, fill_c, outline_c, row, idx)
+                rows.append(row)
+
+            elif geom.geom_type == "MultiPolygon":
+                mp = _make_patch_from_multipolygon(geom)
+                for _, m in enumerate(mp):
+                    mp_copy = row.to_dict()
+                    mp_copy["geometry"] = m
+                    assign_fill_and_outline_to_row(shapes, fill_c, outline_c, mp_copy, idx)
+                    rows.append(mp_copy)
+
+            elif geom.geom_type == "Point":
+                row = row.to_dict()
+                row["geometry"] = Circle((geom.x, geom.y), radius=row["radius"])
+                assign_fill_and_outline_to_row(shapes, fill_c, outline_c, row, idx)
+                rows.append(row)
+
+        patches = pd.DataFrame(rows)
 
         return PatchCollection(
-            patches,
+            patches["geometry"].values.tolist(),
             snap=False,
-            # zorder=4,
             lw=render_params.outline_params.linewidth,
-            facecolor=fill_c,
-            edgecolor=outline_c,
+            facecolor=patches["fill_c"],
+            edgecolor=None if all(outline is None for outline in outline_c) else outline_c,
             **kwargs,
         )
 
@@ -167,6 +217,8 @@ def _get_collection_shape(
     if len(color_vector) == 0:
         color_vector = [render_params.cmap_params.na_color]
 
+    shapes = pd.concat(shapes, ignore_index=True)
+    shapes = gpd.GeoDataFrame(shapes, geometry="geometry")
     _cax = _get_collection_shape(
         shapes=shapes,
         s=render_params.size,
@@ -178,6 +230,7 @@ def _get_collection_shape(
         outline_alpha=render_params.outline_alpha
         # **kwargs,
     )
+
     cax = ax.add_collection(_cax)
 
     # Using dict.fromkeys here since set returns in arbitrary order

src/spatialdata_plot/pl/utils.py

@@ -10,10 +10,13 @@
 from typing import Any, Literal
 
 import matplotlib
+import matplotlib.patches as mpatches
+import matplotlib.path as mpath
 import matplotlib.pyplot as plt
 import multiscale_spatial_image as msi
 import numpy as np
 import pandas as pd
+import shapely
 import spatial_image
 import spatialdata as sd
 import xarray as xr
@@ -299,7 +302,7 @@ def _get_extent_after_transformations(element: Any, cs_name: str) -> Sequence[in
                     if shapes_key == e_id:
 
                         def get_point_bb(
-                            point: Point, radius: int, method: Literal["topleft", "bottomright"], buffer: int = 1
+                            point: Point, radius: int, method: Literal["topleft", "bottomright"], buffer: int = 0
                         ) -> Point:
                             x, y = point.coords[0]
                             if method == "topleft":
@@ -346,7 +349,12 @@ def get_point_bb(
                         del tmp_points
                         del tmp_polygons
 
-                        extent[cs_name][e_id] = x_dims + y_dims
+                        xmin = np.min(x_dims)
+                        xmax = np.max(x_dims)
+                        ymin = np.min(y_dims)
+                        ymax = np.max(y_dims)
+
+                        extent[cs_name][e_id] = [xmin, xmax, ymin, ymax]
 
                         transformations = get_transformation(sdata.shapes[e_id], to_coordinate_system=cs_name)
                         transformations = _flatten_transformation_sequence(transformations)
@@ -1166,6 +1174,51 @@ def _robust_transform(element: Any, cs: str) -> Any:
     return element
 
 
+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 == "Polygon":
+            exterior_coords = geom.exterior.coords[:]
+            interior_coords = []
+            for interior in geom.interiors:
+                interior_coords += interior.coords[:]
+        elif 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"]
+        else:
+            raise ValueError("Unhandled geometry type: " + repr(geom.type))
+
+    return interior_coords, exterior_coords
+
+
+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))]
+
+    return patches
+
+
 def _mpl_ax_contains_elements(ax: Axes) -> bool:
     """Check if any objects have been plotted on the axes object.
 

tests/pl/test_render_shapes.py

@@ -1,8 +1,12 @@
+import anndata
 import geopandas as gpd
 import matplotlib
+import pandas as pd
 import scanpy as sc
 import spatialdata_plot  # noqa: F401
+from shapely.geometry import MultiPolygon, Point, Polygon
 from spatialdata import SpatialData
+from spatialdata.models import ShapesModel, TableModel
 
 from tests.conftest import PlotTester, PlotTesterMeta
 
@@ -57,3 +61,31 @@ def test_plot_can_render_circles_with_default_outline_width(self, sdata_blobs: S
 
     def test_plot_can_render_circles_with_specified_outline_width(self, sdata_blobs: SpatialData):
         sdata_blobs.pl.render_shapes(elements="blobs_circles", outline=True, outline_width=3.0).pl.show()
+
+    def test_plot_can_render_multipolygons(self):
+        def _make_multi():
+            hole = MultiPolygon(
+                [(((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)), [((0.2, 0.2), (0.2, 0.8), (0.8, 0.8), (0.8, 0.2))])]
+            )
+            overlap = MultiPolygon(
+                [
+                    Polygon([(2.0, 0.0), (2.0, 0.8), (2.8, 0.8), (2.8, 0.0)]),
+                    Polygon([(2.2, 0.2), (2.2, 1.0), (3.0, 1.0), (3.0, 0.2)]),
+                ]
+            )
+            poly = Polygon([(4.0, 0.0), (4.0, 1.0), (5.0, 1.0), (5.0, 0.0)])
+            circ = Point(6.0, 0.5)
+            polygon_series = gpd.GeoSeries([hole, overlap, poly, circ])
+            cell_polygon_table = gpd.GeoDataFrame(geometry=polygon_series)
+            sd_polygons = ShapesModel.parse(cell_polygon_table)
+            sd_polygons.loc[:, "radius"] = [None, None, None, 0.3]
+
+            return sd_polygons
+
+        sdata = SpatialData(shapes={"p": _make_multi()})
+        adata = anndata.AnnData(pd.DataFrame({"p": ["hole", "overlap", "square", "circle"]}))
+        adata.obs.loc[:, "region"] = "p"
+        adata.obs.loc[:, "val"] = [1, 2, 3, 4]
+        table = TableModel.parse(adata, region="p", region_key="region", instance_key="val")
+        sdata.table = table
+        sdata.pl.render_shapes(color="val", outline=True, fill_alpha=0.3).pl.show()