Fix axes mismatch when elements have multi-CS transformations (#176)

When elements have transformations to multiple coordinate systems,
filter_by_coordinate_system cannot strip the extra transformations
(upstream limitation). This caused show() to auto-detect too many
CS and raise a mismatch error when the user passed their own axes.

Apply a stricter CS filter when coordinate_systems is auto-detected
and axes are provided: keep only CS that have element types for all
render commands. Also improve the error message with a hint to pass
coordinate_systems= explicitly when the mismatch cannot be resolved.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: timtreis

src/spatialdata_plot/pl/basic.py

@@ -993,9 +993,11 @@ def show(
             ax_x_min, ax_x_max = ax.get_xlim()
             ax_y_max, ax_y_min = ax.get_ylim()  # (0, 0) is top-left
 
-        coordinate_systems = sdata.coordinate_systems if coordinate_systems is None else coordinate_systems
+        cs_was_auto = coordinate_systems is None
+        coordinate_systems = list(sdata.coordinate_systems) if cs_was_auto else coordinate_systems
         if isinstance(coordinate_systems, str):
             coordinate_systems = [coordinate_systems]
+        assert coordinate_systems is not None  # guaranteed by the branches above
 
         for cs in coordinate_systems:
             if cs not in sdata.coordinate_systems:
@@ -1019,14 +1021,43 @@ def show(
             elements=elements_to_be_rendered,
         )
 
+        # When coordinate_systems was auto-detected and the user provided axes,
+        # apply a stricter filter: keep only CS that have element types for ALL
+        # render commands. This handles the case where filter_by_coordinate_system
+        # leaves extra CS due to elements having multi-CS transformations (#176).
+        if ax is not None and cs_was_auto:
+            n_ax = 1 if isinstance(ax, Axes) else len(ax)
+            if len(coordinate_systems) > n_ax:
+                render_type_flags = {
+                    "render_images": "has_images",
+                    "render_labels": "has_labels",
+                    "render_points": "has_points",
+                    "render_shapes": "has_shapes",
+                }
+                required_flags = [render_type_flags[cmd] for cmd in cmds if cmd in render_type_flags]
+                strict_cs = []
+                for cs_name in coordinate_systems:
+                    row = cs_contents.query(f"cs == '{cs_name}'").iloc[0]
+                    if all(row[flag] for flag in required_flags):
+                        strict_cs.append(cs_name)
+                if strict_cs:
+                    coordinate_systems = strict_cs
+
         # catch error in ruff-friendly way
         if ax is not None:  # we'll generate matching number then
             n_ax = 1 if isinstance(ax, Axes) else len(ax)
             if len(coordinate_systems) != n_ax:
-                raise ValueError(
+                msg = (
                     f"Mismatch between number of matplotlib axes objects ({n_ax}) "
                     f"and number of coordinate systems ({len(coordinate_systems)})."
                 )
+                if cs_was_auto:
+                    msg += (
+                        " This can happen when elements have transformations to multiple "
+                        "coordinate systems (e.g. after filter_by_coordinate_system). "
+                        "Pass `coordinate_systems=` explicitly to select which ones to plot."
+                    )
+                raise ValueError(msg)
 
         # set up canvas
         fig_params, scalebar_params = _prepare_params_plot(

src/spatialdata_plot/pl/utils.py

@@ -264,6 +264,7 @@ def _prepare_params_plot(
         if ax is not None and len(ax) != num_panels:
             raise ValueError(f"Len of `ax`: {len(ax)} is not equal to number of panels: {num_panels}.")
         if fig is None:
+            # TODO(#579): infer fig from ax[0].get_figure() instead of requiring it
             raise ValueError(
                 f"Invalid value of `fig`: {fig}. If a list of `Axes` is passed, a `Figure` must also be specified."
             )

tests/pl/test_render.py

@@ -1,5 +1,8 @@
 import matplotlib.pyplot as plt
+import numpy as np
 import pytest
+import spatialdata as sd
+from spatialdata.transformations import Identity
 
 
 def test_render_images_can_plot_one_cyx_image(request):
@@ -62,3 +65,66 @@ def test_keyerror_when_shape_element_does_not_exist(request):
 
     with pytest.raises(KeyError):
         sdata.pl.render_shapes(element="not_found").pl.show()
+
+
+# Regression tests for #176: plotting with user-supplied ax when elements
+# have transformations to multiple coordinate systems.
+
+
+def _make_multi_cs_sdata():
+    """Create sdata where shapes has transformations to two CS but image has only one."""
+    import geopandas as gpd
+    from shapely import Point
+
+    image = sd.models.Image2DModel.parse(
+        np.zeros((1, 10, 10)), dims=("c", "y", "x"), transformations={"aligned": Identity()}
+    )
+    shapes = sd.models.ShapesModel.parse(
+        gpd.GeoDataFrame(geometry=[Point(5, 5)], data={"radius": [2]}),
+        transformations={"aligned": Identity(), "global": Identity()},
+    )
+    return sd.SpatialData(images={"img": image}, shapes={"shp": shapes})
+
+
+def test_single_ax_after_filter_by_coordinate_system():
+    """After filter_by_coordinate_system, single ax should work without specifying CS."""
+    sdata = _make_multi_cs_sdata()
+    sdata_filt = sdata.filter_by_coordinate_system("aligned")
+
+    _, ax = plt.subplots(1, 1)
+    sdata_filt.pl.render_images("img").pl.render_shapes("shp").pl.show(ax=ax)
+    assert ax.get_title() == "aligned"
+
+
+def test_single_ax_with_explicit_cs():
+    """Explicit coordinate_systems with single ax should work."""
+    sdata = _make_multi_cs_sdata()
+
+    _, ax = plt.subplots(1, 1)
+    sdata.pl.render_images("img").pl.render_shapes("shp").pl.show(ax=ax, coordinate_systems="aligned")
+    assert ax.get_title() == "aligned"
+
+
+def test_single_ax_explicit_multi_cs_raises():
+    """Explicitly requesting more CS than axes should still raise."""
+    sdata = _make_multi_cs_sdata()
+
+    _, ax = plt.subplots(1, 1)
+    with pytest.raises(ValueError, match="Mismatch"):
+        sdata.pl.render_shapes("shp").pl.show(ax=ax, coordinate_systems=["aligned", "global"])
+
+
+def test_single_ax_auto_cs_unresolvable_raises():
+    """When strict filtering can't resolve the mismatch, error includes hint."""
+    import geopandas as gpd
+    from shapely import Point
+
+    shapes = sd.models.ShapesModel.parse(
+        gpd.GeoDataFrame(geometry=[Point(5, 5)], data={"radius": [2]}),
+        transformations={"aligned": Identity(), "global": Identity()},
+    )
+    sdata = sd.SpatialData(shapes={"shp": shapes})
+
+    _, ax = plt.subplots(1, 1)
+    with pytest.raises(ValueError, match="coordinate_systems="):
+        sdata.pl.render_shapes("shp").pl.show(ax=ax)