render_shapes cannot handle shapes with multiple holes #505
Description
I encountered an issue when trying to run render_shapes, which said ValueError: 'vertices' must be 2D with shape (N, 2), but your input has shape (10,).
Upon further investigation, I believe that the issue is that one of my shapes has two holes, which the method is not equipped for. Below I am providing a minimal example, the methods are taken directly from spatialdata_plot/pl/utils.py.
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.patches as mpatches
import matplotlib.path as mpath
from shapely import wkt
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
# FROM SPATIALDATA-PLOT
def _split_multipolygon_into_outer_and_inner(mp: 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[:]
else:
raise ValueError(f"Unhandled geometry type: {repr(geom.type)}")
return interior_coords, exterior_coords
def _make_patch_from_multipolygon(mp: 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)
print(f"{vertices=}")
print(f"{all_codes=}")
print(f"{inside=}")
print(f"{outside=}")
patches += [mpatches.PathPatch(mpath.Path(vertices, all_codes))]
return patches
# HELPER METHOD FOR VISUALIZATION
def plot_multipoly(geom, ax=None):
"""
Plot a Polygon or MultiPolygon with automatic coloring for each part.
Each polygon component (if MultiPolygon) is given a distinct color.
Holes are shown as white.
"""
if ax is None:
fig, ax = plt.subplots()
# color cycle from Matplotlib
color_cycle = plt.rcParams["axes.prop_cycle"].by_key()["color"]
def add_patch(geom, ax, color_idx=0):
if isinstance(geom, Polygon):
color = color_cycle[color_idx % len(color_cycle)]
patch = mpatches.Polygon(
np.array(geom.exterior.coords),
closed=True,
facecolor=color,
edgecolor="black",
alpha=0.5,
)
ax.add_patch(patch)
# draw holes as white patches
for interior in geom.interiors:
hole_patch = mpatches.Polygon(
np.array(interior.coords),
closed=True,
facecolor="white",
edgecolor="black",
)
ax.add_patch(hole_patch)
elif isinstance(geom, MultiPolygon):
for i, poly in enumerate(geom.geoms):
add_patch(poly, ax, color_idx=i)
else:
raise ValueError(f"Unsupported geometry type: {type(geom)}")
add_patch(geom, ax)
# adjust plot limits
minx, miny, maxx, maxy = geom.bounds
ax.set_xlim(minx - 1, maxx + 1)
ax.set_ylim(miny - 1, maxy + 1)
ax.set_aspect("equal")
plt.show()
# the polygon that originaly triggered the issue
# geom = wkt.loads("MULTIPOLYGON (((545.5 78.5, 543.5 78.5, 543.5 79.5, 542.5 79.5, 542.5 81.5, 544.5 81.5, 543.5 85.5, 544.5 85.5, 544.5 88.5, 545.5 88.5, 545.5 89.5, 549.5 90.5, 549.5 89.5, 547.5 88.5, 547.5 84.5, 545.5 84.5, 545.5 81.5, 550.5 81.5, 550.5 77.5, 549.5 77.5, 549.5 76.5, 548.5 76.5, 548.5 77.5, 547.5 77.5, 547.5 76.5, 545.5 76.5, 545.5 78.5), (547.5 78.5, 546.5 78.5, 546.5 77.5, 547.5 77.5, 547.5 78.5), (549.5 78.5, 549.5 79.5, 547.5 79.5, 547.5 78.5, 549.5 78.5)), ((550.5 89.5, 550.5 88.5, 549.5 88.5, 549.5 89.5, 550.5 89.5)))")
geom = wkt.loads("MULTIPOLYGON (((0 0, 5 0, 5 5, 0 5, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1), (3 3, 3 4, 4 4, 4 3, 3 3)))")
# print the polygon string
print(geom)
# visualizing the polygon
plot_multipoly(geom)
# using the methods from spatialdata-plot
_make_patch_from_multipolygon(geom)
When rendered properly, the polygon should look like this:
I get the following output using your methods:
vertices=array([[0., 0.],
[5., 0.],
[5., 5.],
[0., 5.],
[0., 0.],
[3., 3.],
[4., 3.],
[4., 4.],
[3., 4.],
[3., 3.],
[1., 1.],
[1., 2.],
[2., 2.],
[2., 1.],
[1., 1.]])
all_codes=array([1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2],
dtype=uint8)
inside=[(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0), (3.0, 3.0), (3.0, 4.0), (4.0, 4.0), (4.0, 3.0), (3.0, 3.0)]
outside=[(0.0, 0.0), (5.0, 0.0), (5.0, 5.0), (0.0, 5.0), (0.0, 0.0)]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[65], line 119
117 plot_multipoly(geom)
118 # using the methods from spatialdata-plot
--> 119 _make_patch_from_multipolygon(geom)
Cell In[65], line 55, in _make_patch_from_multipolygon(mp)
52 print(f"{inside=}")
53 print(f"{outside=}")
---> 55 patches += [mpatches.PathPatch(mpath.Path(vertices, all_codes))]
57 return patches
File /g/huber/users/meyerben/notebooks/spatial_transcriptomics/SegTraQ/.venv/lib/python3.13/site-packages/matplotlib/path.py:135, in Path.__init__(self, vertices, codes, _interpolation_steps, closed, readonly)
133 codes = np.asarray(codes, self.code_type)
134 if codes.ndim != 1 or len(codes) != len(vertices):
--> 135 raise ValueError("'codes' must be a 1D list or array with the "
136 "same length of 'vertices'. "
137 f"Your vertices have shape {vertices.shape} "
138 f"but your codes have shape {codes.shape}")
139 if len(codes) and codes[0] != self.MOVETO:
140 raise ValueError("The first element of 'code' must be equal "
141 f"to 'MOVETO' ({self.MOVETO}). "
142 f"Your first code is {codes[0]}")
ValueError: 'codes' must be a 1D list or array with the same length of 'vertices'. Your vertices have shape (15, 2) but your codes have shape (20,)
I have also provided the second polygon in the code (that triggered the original issue).