Document that spatial queries on MultiPolygon keep also the components outside the queried region

[LucaMarconato]

Bug report from @claude, bug discovered and triaged by @LucaMarconato.

Spatial queries on MultiPolygon shapes return the whole geometry, including the components that lie entirely outside the queried region.

MWE

blobs()["blobs_multipolygons"] holds 2 MultiPolygons, each made of 2 disjoint components. Querying a region that only contains the first component of each returns both MultiPolygons unchanged:

# /// script
# requires-python = ">=3.12"
# dependencies = [
#     "spatialdata @ git+https://github.com/scverse/spatialdata.git@65dc73e",
#     "spatialdata-plot",
#     "matplotlib",
# ]
# ///
"""Spatial queries on MultiPolygon keep also the components that lie fully outside the queried region.

`blobs()["blobs_multipolygons"]` contains 2 MultiPolygons, each made of 2 disjoint Polygon components.
Querying a region that only contains the *first* component of each MultiPolygon returns the MultiPolygons
unchanged: the second component - which has zero overlap with the query region - is still there.

`polygon_query(..., clip=True)` does drop them, but `bounding_box_query()` has no `clip` argument at all,
so there is no way to get the same result from a bounding-box query.
"""

import warnings

import matplotlib.pyplot as plt
import spatialdata_plot  # noqa: F401  (registers the .pl accessor)
from matplotlib.patches import Rectangle
from shapely.geometry import box

from spatialdata import SpatialData, bounding_box_query, polygon_query
from spatialdata.datasets import blobs

warnings.simplefilter("ignore")

# ---------------------------------------------------------------- data + query
sdata = blobs()
multipolygons = sdata["blobs_multipolygons"]
XMIN, YMIN, XMAX, YMAX = 285.0, 190.0, 345.0, 262.0
query_region = box(XMIN, YMIN, XMAX, YMAX)

no_clip = polygon_query(multipolygons, polygon=query_region, target_coordinate_system="global", clip=False)
clipped = polygon_query(multipolygons, polygon=query_region, target_coordinate_system="global", clip=True)
bbox = bounding_box_query(
    multipolygons,
    axes=("x", "y"),
    min_coordinate=[XMIN, YMIN],
    max_coordinate=[XMAX, YMAX],
    target_coordinate_system="global",
)  # no `clip` argument available here


def describe(name, gdf):
    print(f"{name}:")
    n_outside = 0
    for idx, geom in zip(gdf.index, gdf.geometry, strict=True):
        parts = list(geom.geoms) if geom.geom_type == "MultiPolygon" else [geom]
        outside = [p for p in parts if not p.intersects(query_region)]
        n_outside += len(outside)
        print(
            f"  index {idx}: {geom.geom_type} with {len(parts)} component(s), "
            f"{len(outside)} of which do NOT intersect the query region"
        )
    return n_outside


print(f"query region: box({XMIN}, {YMIN}, {XMAX}, {YMAX})\n")
describe("original", multipolygons)
n_outside = describe("polygon_query(..., clip=False)  <-- default", no_clip)
describe("polygon_query(..., clip=True)", clipped)
n_outside += describe("bounding_box_query(...)         <-- no `clip` argument exists", bbox)
print(
    "\nexpected: the MultiPolygon components that do not intersect the query region are dropped,\n"
    "or at least this behaviour (and the `clip=True` workaround) is documented\n"
    f"VERDICT: {'BUG REPRODUCED' if n_outside else 'NOT REPRODUCED'}"
)

# ------------------------------------------------------------------- plotting
panels = {
    "original": multipolygons,
    "polygon_query(clip=False), default\nand bounding_box_query():\ncomponents outside the box are kept": no_clip,
    "polygon_query(clip=True):\ncomponents outside the box are dropped": clipped,
}
fig, axes = plt.subplots(1, 3, figsize=(15, 5.5))
for ax, (title, shapes) in zip(axes, panels.items(), strict=True):
    SpatialData(images={"blobs_image": sdata["blobs_image"]}, shapes={"shapes": shapes}).pl.render_images(
        "blobs_image"
    ).pl.render_shapes("shapes", fill_alpha=0.6, outline_alpha=1.0, outline_color="white").pl.show(ax=ax, title=title)
    ax.add_patch(
        Rectangle((XMIN, YMIN), XMAX - XMIN, YMAX - YMIN, fill=False, edgecolor="red", lw=2, ls="--", zorder=10)
    )
fig.tight_layout()
fig.savefig("multipolygon_query_keeps_outside_components.png", dpi=120)
print("figure written to multipolygon_query_keeps_outside_components.png")

Full self-contained repro (PEP 723, uv run repro.py) attached; it also renders the figure below with spatialdata-plot.

The left and middle panels are identical: the default query returns everything. Red dashed box = query region.

Why this is confusing

For a single Polygon that only partially overlaps the region, returning the whole geometry is the documented, intended behaviour ("keep the shape if it intersects"). For a MultiPolygon the same rule silently keeps components that have zero overlap with the query region, which is much more surprising: the returned element can extend arbitrarily far outside the queried region.

clip=True gives the expected result, but:

• it is not obvious that this is the knob to reach for — the parameter reads as "trim the boundary", not "drop the parts that are not in the region";
• it is only available in polygon_query; bounding_box_query has no clip argument, so there is no way to get the same result from a bounding-box query;
• as a side effect it changes the geometry type (MultiPolygon -> Polygon) when only one component survives.

Suggestion

Mainly a documentation issue: document explicitly, in the polygon_query/bounding_box_query docstrings and in the query section of the docs, that shapes are returned whole and that for MultiPolygon this includes components lying fully outside the queried region, and point to clip=True as the way to get geometrically cropped output.

Additionally, consider adding clip to bounding_box_query (it could simply forward to the polygon-query path with a box), so that the workaround is available for both query types.

Environment

uv run repro.py with the PEP 723 metadata in the script (fresh, isolated environment; spatialdata built from main @ 65dc73e; Python 3.13, latest releases of the dependencies at run time). macOS (arm64).

---

Issue generated by Claude.