diff --git a/src/spatialdata_plot/pl/_validate.py b/src/spatialdata_plot/pl/_validate.py
index a979ad02..b538bfc4 100644
--- a/src/spatialdata_plot/pl/_validate.py
+++ b/src/spatialdata_plot/pl/_validate.py
@@ -99,6 +99,7 @@ def _validate_show_parameters(
     xlabel: str | None,
     ylabel: str | None,
     pad_extent: int | float,
+    crop_coord: tuple[float, float, float, float] | None,
     ax: list[Axes] | Axes | None,
     return_ax: bool,
     save: str | Path | None,
@@ -192,6 +193,20 @@ def _validate_show_parameters(
     if not isinstance(pad_extent, int | float):
         raise TypeError("Parameter 'pad_extent' must be numeric.")
 
+    if crop_coord is not None:
+        if (
+            not isinstance(crop_coord, tuple)
+            or len(crop_coord) != 4
+            or not all(isinstance(v, int | float) for v in crop_coord)
+        ):
+            raise TypeError("Parameter 'crop_coord' must be a tuple of four numbers (xmin, xmax, ymin, ymax).")
+        xmin, xmax, ymin, ymax = crop_coord
+        if not (xmin < xmax and ymin < ymax):
+            raise ValueError(
+                f"Parameter 'crop_coord' must satisfy xmin < xmax and ymin < ymax, got (xmin={xmin}, xmax={xmax}, "
+                f"ymin={ymin}, ymax={ymax})."
+            )
+
     if ax is not None and not isinstance(ax, Axes | list):
         raise TypeError("Parameter 'ax' must be a matplotlib.axes.Axes or a list of Axes.")
 
diff --git a/src/spatialdata_plot/pl/basic.py b/src/spatialdata_plot/pl/basic.py
index b910671c..4eaab105 100644
--- a/src/spatialdata_plot/pl/basic.py
+++ b/src/spatialdata_plot/pl/basic.py
@@ -56,6 +56,7 @@
     _split_colorbar_params,
 )
 from spatialdata_plot.pl.render_params import (
+    BBox,
     CBAR_DEFAULT_FRACTION,
     CBAR_DEFAULT_LOCATION,
     CBAR_DEFAULT_PAD,
@@ -1318,6 +1319,7 @@ def show(
         xlabel: str | None = None,
         ylabel: str | None = None,
         pad_extent: int | float = 0,
+        crop_coord: tuple[float, float, float, float] | None = None,
         ax: list[Axes] | Axes | None = None,
         return_ax: bool = False,
         save: str | Path | None = None,
@@ -1385,7 +1387,21 @@ def show(
         ylabel : str | None, default None
             Label for the y axis, applied to every rendered panel. ``None`` leaves it unlabelled.
         pad_extent : int | float, default 0
-            Padding added around the computed spatial extent on all sides.
+            Padding added around the computed spatial extent on all sides. Ignored when ``crop_coord`` is set.
+        crop_coord : tuple[float, float, float, float] | None, default None
+            Restrict the plot to a bounding box ``(xmin, xmax, ymin, ymax)`` in the rendered coordinate
+            system's units (same order as :meth:`matplotlib.axes.Axes.axis` and as ``crop_coord`` in
+            :func:`scanpy.pl.spatial` / :func:`squidpy.pl.spatial_scatter`). Points and shapes are
+            subsetted before drawing for speed; large images (both the matplotlib and ``datashader``
+            backends) and labels are rasterized to the window only, so the full array is never
+            materialized (fast at Visium HD scale) and the zoom keeps full figure resolution. For points
+            and shapes, auto-scaled color ranges are computed from the full element so colors match the
+            uncropped plot; a cropped image's contrast auto-scales over the window (pass explicit
+            ``vmin``/``vmax`` or a ``norm`` to fix it). Labels coloured by a plain-string column with the
+            default palette fall back to full render + clip (windowing could otherwise reshuffle
+            colours); use a categorical dtype or an explicit palette to keep the windowed fast path.
+            Requires a single coordinate system (pass ``coordinate_systems`` with one entry if several
+            would otherwise be rendered).
         ax : list[Axes] | Axes | None
             Pre-existing matplotlib axes to plot on. Can be a single :class:`~matplotlib.axes.Axes` or a list
             matching the number of coordinate systems. If ``None``, a new figure and axes are created.
@@ -1449,6 +1465,7 @@ def show(
             xlabel=xlabel,
             ylabel=ylabel,
             pad_extent=pad_extent,
+            crop_coord=crop_coord,
             ax=ax,
             return_ax=return_ax,
             save=save,
@@ -1502,6 +1519,19 @@ def show(
             ax=ax,
         )
 
+        # `crop_coord` is one box in one coordinate system's units; applying the same numbers across
+        # coordinate systems with different scales/units would be incoherent. Convert the public
+        # (xmin, xmax, ymin, ymax) order into the internal BBox once, here at the boundary.
+        crop_box: BBox | None = None
+        if crop_coord is not None:
+            if len(coordinate_systems) > 1:
+                raise ValueError(
+                    f"`crop_coord` requires a single coordinate system, but {len(coordinate_systems)} would be "
+                    f"rendered ({coordinate_systems}). Pass `coordinate_systems=` with exactly one entry."
+                )
+            xmin, xmax, ymin, ymax = crop_coord
+            crop_box = BBox(xmin, ymin, xmax, ymax)
+
         panels = _plan_panels(
             coordinate_systems=coordinate_systems,
             render_cmds=render_cmds,
@@ -1534,6 +1564,7 @@ def show(
             scalebar_units=scalebar_units,
             scalebar_kwargs=scalebar_params,
         )
+        fig_params.crop = crop_box
         legend_params_obj = _build_legend_params(
             legend_params=legend_params,
             legend_fontsize=legend_fontsize,
@@ -1591,25 +1622,33 @@ def show(
                         "all geometries are empty. Drop the element or restore at least one non-empty geometry."
                     )
 
-            # fast path for axis-aligned transforms; identical result, falls back to get_extent otherwise
-            extent = _get_extent_fast(
-                sdata,
-                coordinate_system=cs,
-                has_images=has_images and wants["images"],
-                has_labels=has_labels and wants["labels"],
-                has_points=has_points and wants["points"],
-                has_shapes=has_shapes and wants["shapes"],
-                elements=wanted_elements,
-            )
-            cs_x_min, cs_x_max = extent["x"]
-            cs_y_min, cs_y_max = extent["y"]
-
-            if any([has_images, has_labels, has_points, has_shapes]):
-                # If the axis already has limits, only expand them but not overwrite
-                x_min = min(ax_x_min, cs_x_min) - pad_extent
-                x_max = max(ax_x_max, cs_x_max) + pad_extent
-                y_min = min(ax_y_min, cs_y_min) - pad_extent
-                y_max = max(ax_y_max, cs_y_max) + pad_extent
+            if crop_box is not None:
+                # `crop_coord` pins the view to the exact box: bypass the expand-don't-overwrite merge
+                # (which would expand back out to any pre-existing axes limits) and ignore `pad_extent`.
+                x_min, x_max, y_min, y_max = crop_box.x0, crop_box.x1, crop_box.y0, crop_box.y1
+                set_limits = True
+            else:
+                # fast path for axis-aligned transforms; identical result, falls back to get_extent otherwise
+                extent = _get_extent_fast(
+                    sdata,
+                    coordinate_system=cs,
+                    has_images=has_images and wants["images"],
+                    has_labels=has_labels and wants["labels"],
+                    has_points=has_points and wants["points"],
+                    has_shapes=has_shapes and wants["shapes"],
+                    elements=wanted_elements,
+                )
+                cs_x_min, cs_x_max = extent["x"]
+                cs_y_min, cs_y_max = extent["y"]
+                set_limits = any([has_images, has_labels, has_points, has_shapes])
+                if set_limits:
+                    # If the axis already has limits, only expand them but not overwrite
+                    x_min = min(ax_x_min, cs_x_min) - pad_extent
+                    x_max = max(ax_x_max, cs_x_max) + pad_extent
+                    y_min = min(ax_y_min, cs_y_min) - pad_extent
+                    y_max = max(ax_y_max, cs_y_max) + pad_extent
+
+            if set_limits:
                 ax.set_xlim(x_min, x_max)
                 ax.set_ylim(y_max, y_min)  # (0, 0) is top-left
 
diff --git a/src/spatialdata_plot/pl/render.py b/src/spatialdata_plot/pl/render.py
index d1242559..b2aa0a51 100644
--- a/src/spatialdata_plot/pl/render.py
+++ b/src/spatialdata_plot/pl/render.py
@@ -2,7 +2,7 @@
 
 import dataclasses
 from collections import abc
-from collections.abc import Sequence
+from collections.abc import Callable, Sequence
 from copy import copy
 from typing import Any, Literal, cast
 
@@ -78,6 +78,9 @@
     colormap_with_alpha,
 )
 from spatialdata_plot.pl.utils import (
+    _bbox_mask_points,
+    _bbox_mask_shapes,
+    _bbox_to_element_space,
     _categorical_legend_handles,
     _decorate_axs,
     _fast_extent,
@@ -105,6 +108,55 @@
 )
 
 
+def _full_element_norm_range(
+    color_spec: ColorSpec, cmap_params: CmapParams, transfunc: Callable[..., Any] | None
+) -> tuple[float, float] | None:
+    """``(vmin, vmax)`` over the full (pre-crop) continuous color vector, after ``transfunc``.
+
+    A crop drops rows before drawing, but the matplotlib backend must keep the uncropped colour range
+    so the zoom matches the full plot. Returns the range resolved over the full element's transfunc'd
+    values (what the uncropped plot autoscales to), or ``None`` for categorical color or when the user
+    already fixed vmin/vmax. Non-mutating and applied only on the matplotlib path: the datashader and
+    image backends never receive it and autoscale over the visible window instead.
+    """
+    if not color_spec.is_continuous:
+        return None
+    if cmap_params.norm.vmin is not None and cmap_params.norm.vmax is not None:
+        return (cmap_params.norm.vmin, cmap_params.norm.vmax)
+    vec = color_spec.apply_transfunc(transfunc).color_vector if transfunc is not None else color_spec.color_vector
+    resolved = _resolve_continuous_norm(vec, cmap_params)
+    return (float(resolved.vmin), float(resolved.vmax))
+
+
+def _pin_norm_to_full_range(
+    norm: Normalize,
+    full_color_spec: ColorSpec | None,
+    cmap_params: CmapParams,
+    transfunc: Callable[..., Any] | None,
+) -> None:
+    """Pin ``norm`` in place to the full (pre-crop) continuous range so a cropped matplotlib zoom keeps
+    the uncropped plot's colours. No-op when uncropped (``full_color_spec is None``) or categorical.
+    """
+    if full_color_spec is None:
+        return
+    full_range = _full_element_norm_range(full_color_spec, cmap_params, transfunc)
+    if full_range is not None:
+        norm.vmin, norm.vmax = full_range
+
+
+def _crop_color_is_stable(color_col: pd.Series | None, palette: Any) -> bool:
+    """Whether windowing a labels raster keeps colours stable (else the caller full-renders + clips).
+
+    Windowing drops off-window instances before colours resolve. That is colour-stable for numeric
+    (continuous) columns, pandas-Categorical columns (fixed levels), and explicit palettes; a plain-string
+    column with the default palette would reshuffle colours as whole categories drop out of the window
+    (their sorted positions shift), so windowing is unsafe there.
+    """
+    if palette is not None or color_col is None:
+        return True
+    return pd.api.types.is_numeric_dtype(color_col) or isinstance(color_col.dtype, pd.CategoricalDtype)
+
+
 def _get_top_data_array(element: xr.DataArray | DataTree) -> xr.DataArray:
     if isinstance(element, DataTree):
         return next(iter(next(iter(element.values())).data_vars.values()))
@@ -700,6 +752,23 @@ def _render_shapes(
         if outline_color_spec is not None:
             outline_color_spec = outline_color_spec.filter(keep)
 
+    # On-the-fly crop: drop shapes outside the box. Keep a reference to the pre-crop colour spec so the
+    # matplotlib backend can pin the uncropped norm; the range itself is resolved lazily below (only on
+    # the matplotlib path — datashader autoscales over the window, so it must not pay the full scan).
+    full_color_spec: ColorSpec | None = None
+    if fig_params.crop is not None:
+        elem_bbox = _bbox_to_element_space(sdata_filt.shapes[element], coordinate_system, fig_params.crop)
+        if elem_bbox is not None:  # None = rotation/shear; render full and let axis limits clip
+            full_color_spec = color_spec  # capture before filtering (filter returns a new object)
+            keep_crop = _bbox_mask_shapes(shapes, elem_bbox)
+            color_spec = color_spec.filter(keep_crop)
+            shapes = shapes[keep_crop].reset_index(drop=True)
+            if len(shapes) == 0:
+                return
+            sdata_filt[element] = shapes
+            if outline_color_spec is not None:
+                outline_color_spec = outline_color_spec.filter(keep_crop)
+
     color_spec = color_spec.apply_transfunc(render_params.transfunc)
 
     norm = render_params.cmap_params.fresh_norm()
@@ -917,6 +986,10 @@ def _draw_centroids(xy: np.ndarray, radius: float | None = None) -> None:
         cax = _build_ds_colorbar(reduction_bounds, norm, render_params.cmap_params.cmap)
 
     elif method == "matplotlib":
+        # Under crop, pin the norm to the full-element range so the zoom matches the uncropped plot
+        # (datashader, handled above, autoscales over the window instead). Resolved here (not in the crop
+        # block) so the datashader path never pays the full-element scan.
+        _pin_norm_to_full_range(norm, full_color_spec, render_params.cmap_params, render_params.transfunc)
         # Build the paths once and share them across the fill and outline collections (geometry is
         # identical; only colours/alpha/linewidth differ), then apply the coordinate-system affine
         # once to the shared Path objects rather than once per collection.
@@ -1451,6 +1524,22 @@ def _render_points(
         points = points[keep].reset_index(drop=True)
         _reparse_points(sdata_filt, element, points, transformation_in_cs, coordinate_system, col_for_color)
 
+    # On-the-fly crop: drop points outside the box. Keep a reference to the pre-crop colour spec so the
+    # matplotlib backend can pin the uncropped norm; the range itself is resolved lazily below (only on
+    # the matplotlib path — datashader autoscales over the window, so it must not pay the full scan).
+    full_color_spec: ColorSpec | None = None
+    if fig_params.crop is not None:
+        elem_bbox = _bbox_to_element_space(sdata.points[element], coordinate_system, fig_params.crop)
+        if elem_bbox is not None:  # None = rotation/shear; render full and let axis limits clip
+            keep_crop = _bbox_mask_points(points["x"].to_numpy(), points["y"].to_numpy(), elem_bbox)
+            if not keep_crop.any():
+                return
+            full_color_spec = color_spec  # capture before filtering (filter returns a new object)
+            color_spec = color_spec.filter(keep_crop)
+            points = points[keep_crop].reset_index(drop=True)
+            n_points = len(points)  # method auto-threshold below should see the drawn count
+            _reparse_points(sdata_filt, element, points, transformation_in_cs, coordinate_system, col_for_color)
+
     color_spec = color_spec.apply_transfunc(render_params.transfunc)
 
     trans, trans_data = _prepare_transformation(sdata.points[element], coordinate_system, ax)
@@ -1509,11 +1598,12 @@ def _render_points(
     elif method == "matplotlib":
         # matplotlib colors each point by its own value, so resolve the norm to match shapes/labels
         # instead of letting ax.scatter autoscale a fresh one. Non-continuous keeps the fresh norm.
-        norm = (
-            _resolve_continuous_norm(color_spec.color_vector, render_params.cmap_params)
-            if color_spec.is_continuous
-            else render_params.cmap_params.fresh_norm()
-        )
+        # Under crop, pin to the full-element range so the zoom matches the uncropped plot.
+        if color_spec.is_continuous:
+            norm = _resolve_continuous_norm(color_spec.color_vector, render_params.cmap_params)
+            _pin_norm_to_full_range(norm, full_color_spec, render_params.cmap_params, render_params.transfunc)
+        else:
+            norm = render_params.cmap_params.fresh_norm()
         # update axis limits if plot was empty before (necessary if datashader comes after)
         update_parameters = not _mpl_ax_contains_elements(ax)
         cax = _scatter_points(
@@ -1756,6 +1846,8 @@ def _render_images(
             width=fig_params.fig.get_size_inches()[0],
             height=fig_params.fig.get_size_inches()[1],
             scale=scale,
+            crop=fig_params.crop,
+            extent=extent,
         )
     # rasterize spatial image if necessary to speed up performance
     use_datashader = render_params.method == "datashader"
@@ -1774,6 +1866,7 @@ def _render_images(
             coordinate_system=coordinate_system,
             extent=extent,
             downsample_method=downsample_method,
+            crop=fig_params.crop,
         )
     elif rasterize:
         img = _rasterize_if_necessary(
@@ -1783,6 +1876,7 @@ def _render_images(
             height=fig_params.fig.get_size_inches()[1],
             coordinate_system=coordinate_system,
             extent=extent,
+            crop=fig_params.crop,
         )
 
     channels = img.coords["c"].values.tolist() if render_params.channel is None else render_params.channel
@@ -2182,6 +2276,14 @@ def _render_labels(
     _guard_2d_only(label, element, "labels")
     extent = get_extent(label, coordinate_system=coordinate_system)
 
+    # Only window the label raster when doing so keeps colours stable; otherwise fall back to full
+    # render + axis-clip so cropped colours always match the uncropped plot (see _crop_color_is_stable).
+    crop_labels = fig_params.crop
+    if crop_labels is not None:
+        color_col = sdata[table_name].obs.get(col_for_color) if col_for_color and table_name else None
+        if not _crop_color_is_stable(color_col, palette):
+            crop_labels = None
+
     # get best scale out of multiscale label
     if isinstance(label, DataTree):
         label = _multiscale_to_spatial_image(
@@ -2191,6 +2293,8 @@ def _render_labels(
             height=fig_params.fig.get_size_inches()[1],
             scale=scale,
             is_label=True,
+            crop=crop_labels,
+            extent=extent,
         )
 
     # spatialdata >= 0.8 rejects non-integer label rasters at the model boundary, but the library
@@ -2212,6 +2316,7 @@ def _render_labels(
             height=fig_params.fig.get_size_inches()[1],
             coordinate_system=coordinate_system,
             extent=extent,
+            crop=crop_labels,
         )
 
         # the above adds a useless c dimension of 1 (y, x) -> (1, y, x)
@@ -2221,6 +2326,11 @@ def _render_labels(
     # rasterize mask below; compute them once over the (possibly rasterized) raster.
     unique_labels = np.unique(label.values)
 
+    if crop_labels is not None and not np.any(unique_labels != 0):
+        # The crop window covers no labels (only background). Nothing to draw; the panel's axis
+        # limits already clip to the box, so return before the instance-overlap/colour machinery.
+        return
+
     if table_name is None:
         instance_id = unique_labels
         table = None
diff --git a/src/spatialdata_plot/pl/render_params.py b/src/spatialdata_plot/pl/render_params.py
index f2fa8716..ec09b240 100644
--- a/src/spatialdata_plot/pl/render_params.py
+++ b/src/spatialdata_plot/pl/render_params.py
@@ -3,7 +3,7 @@
 from collections.abc import Callable, Mapping, Sequence
 from copy import copy
 from dataclasses import dataclass, field
-from typing import Any, Literal
+from typing import Any, Literal, NamedTuple
 
 import numpy as np
 from matplotlib.axes import Axes
@@ -16,6 +16,20 @@
 _DsReduction = Literal["sum", "mean", "any", "count", "std", "var", "max", "min"]
 _ImageDsReduction = Literal["max", "min", "mean", "mode", "first", "last", "var", "std"]
 
+
+class BBox(NamedTuple):
+    """Axis-aligned bounding box in one coordinate space, corners sorted so ``x0 <= x1`` and ``y0 <= y1``.
+
+    On-the-fly cropping uses a single ``(x0, y0, x1, y1)`` ordering everywhere internally, so no crop
+    helper has to restate which order it received. ``show()`` converts the public
+    ``crop_coord=(xmin, xmax, ymin, ymax)`` into this once at the boundary.
+    """
+
+    x0: float
+    y0: float
+    x1: float
+    y1: float
+
 # Canonical definition for the package; imported by basic.py and utils.py.
 # replace with
 # from spatialdata._types import ColorLike
@@ -195,6 +209,7 @@ class FigParams:
     title: str | Sequence[str] | None = None
     ax_labels: Sequence[str] | None = None
     frameon: bool | None = None
+    crop: BBox | None = None
 
 
 @dataclass
diff --git a/src/spatialdata_plot/pl/utils.py b/src/spatialdata_plot/pl/utils.py
index 3afff5f4..dc534295 100644
--- a/src/spatialdata_plot/pl/utils.py
+++ b/src/spatialdata_plot/pl/utils.py
@@ -57,6 +57,7 @@
 from spatialdata_plot._logging import logger
 from spatialdata_plot.pl._scanpy_compat import _add_categorical_legend, default_102
 from spatialdata_plot.pl.render_params import (
+    BBox,
     Color,
     ColorbarSpec,
     FigParams,
@@ -805,6 +806,24 @@ def _get_valid_cs(
     return valid_cs
 
 
+def _rasterize_to_bbox(
+    image: DataArray,
+    bbox: BBox,
+    coordinate_system: str,
+    target_x_dims: float,
+    target_y_dims: float,
+) -> DataArray:
+    """Rasterize ``image`` to the world ``bbox`` at ~``target_*_dims`` figure pixels.
+
+    ``rasterize`` interprets ``target_unit_to_pixels`` in world units (not intrinsic pixels), so
+    dividing the target pixel count by the world span keeps placement correct for any transformation
+    (translation, scale, etc.) without a transform rewrite.
+    """
+    x0, y0, x1, y1 = bbox
+    target_unit_to_pixels = min(target_y_dims / (y1 - y0), target_x_dims / (x1 - x0))
+    return rasterize(image, ("y", "x"), [y0, x0], [y1, x1], coordinate_system, target_unit_to_pixels=target_unit_to_pixels)
+
+
 def _rasterize_if_necessary(
     image: DataArray,
     dpi: float,
@@ -812,6 +831,7 @@ def _rasterize_if_necessary(
     height: float,
     coordinate_system: str,
     extent: dict[str, tuple[float, float]],
+    crop: BBox | None = None,
 ) -> DataArray:
     """Ensure fast rendering by adapting the resolution if necessary.
 
@@ -849,35 +869,78 @@ def _rasterize_if_necessary(
     target_y_dims = dpi * height
     target_x_dims = dpi * width
 
-    # Rasterize when the source image is substantially larger than what the
-    # current figure DPI × size requires.  The +100 margin avoids rasterizing
-    # when the image is only slightly larger than the target.
-    do_rasterization = y_dims > target_y_dims + 100 or x_dims > target_x_dims + 100
-
-    if do_rasterization:
-        logger.info("Rasterizing image for faster rendering.")
-        # ``rasterize`` interprets ``target_unit_to_pixels`` in world units, not
-        # intrinsic pixels. Dividing by world extent keeps the result correct
-        # for any transformation (translation, scale, etc.).
-        world_x = float(extent["x"][1]) - float(extent["x"][0])
-        world_y = float(extent["y"][1]) - float(extent["y"][0])
-        target_unit_to_pixels = min(target_y_dims / world_y, target_x_dims / world_x)
-        image = rasterize(
-            image,
-            ("y", "x"),
-            [extent["y"][0], extent["x"][0]],
-            [extent["y"][1], extent["x"][1]],
-            coordinate_system,
-            target_unit_to_pixels=target_unit_to_pixels,
-        )
-        if hasattr(image.data, "compute"):
-            # rasterize is lazy; downstream reads the result once per channel (NaN check,
-            # compositing, draw), so materialize once instead of re-running the warp each time.
-            image = image.copy(data=image.data.compute())
+    if crop is not None:
+        # Restrict rasterization to the crop window: only its source pixels are read (cost bounded by
+        # the window, not the full image), and the zoom keeps full figure resolution.
+        bbox = crop
+    elif y_dims > target_y_dims + 100 or x_dims > target_x_dims + 100:
+        # Rasterize when the source is substantially larger than the figure DPI × size requires; the
+        # +100 margin avoids rasterizing when the image is only slightly larger than the target.
+        bbox = BBox(extent["x"][0], extent["y"][0], extent["x"][1], extent["y"][1])
+    else:
+        return image
 
+    logger.info("Rasterizing image for faster rendering.")
+    image = _rasterize_to_bbox(image, bbox, coordinate_system, target_x_dims, target_y_dims)
+    if hasattr(image.data, "compute"):
+        # rasterize is lazy; downstream reads the result once per channel (NaN check, compositing,
+        # draw), so materialize once instead of re-running the warp each time.
+        image = image.copy(data=image.data.compute())
     return image
 
 
+def _datashader_window_image(
+    image: DataArray,
+    crop: BBox,
+    coordinate_system: str,
+    target_x_dims: int,
+    target_y_dims: int,
+    downsample_method: str,
+) -> DataArray | None:
+    """Datashader-aggregate only the crop window of an image; ``None`` to fall back to the full path.
+
+    ``base`` and the aggregate are made to cover the *same* integer pixel window (the fractional crop
+    box, clamped to the image), so the aggregated grid is index-assigned onto ``base`` without a
+    sub-pixel shift: ``base`` is rasterized over that window mapped back to world (correct coords /
+    transform) and datashader aggregates the sliced source over its own pixel extent. The full image is
+    never read. Returns ``None`` for rotation/shear (box does not map to a box) or an empty window.
+    Images always carry a ``c`` dim (``Image2DModel``) and this path is images-only, so no 2D branch.
+    """
+    y_dims, x_dims = image.shape[1], image.shape[2]
+    elem = _bbox_to_element_space(image, coordinate_system, crop)
+    if elem is None:
+        return None
+    # Integer pixel window: floor/ceil brackets the fractional element box, clamped to the image bounds.
+    px0, py0 = max(int(np.floor(elem.x0)), 0), max(int(np.floor(elem.y0)), 0)
+    px1, py1 = min(int(np.ceil(elem.x1)), x_dims), min(int(np.ceil(elem.y1)), y_dims)
+    if px1 <= px0 or py1 <= py0:
+        return None
+
+    # World box of that integer window (forward affine on its corners), so ``base`` covers exactly the
+    # region the aggregate does — not the fractional crop box, which would shift and (when clamped at an
+    # edge) overrun the sliced source.
+    matrix = get_transformation(image, get_all=True)[coordinate_system].to_affine_matrix(("x", "y"), ("x", "y"))
+    corners = np.array([[px0, py0], [px1, py0], [px0, py1], [px1, py1]], dtype=float)
+    world = corners @ matrix[:2, :2].T + matrix[:2, 2]
+    world_box = BBox(world[:, 0].min(), world[:, 1].min(), world[:, 0].max(), world[:, 1].max())
+
+    base = _rasterize_to_bbox(image, world_box, coordinate_system, target_x_dims, target_y_dims)
+    if hasattr(base.data, "compute"):
+        base = base.copy(data=base.data.compute())  # materialize before we overwrite .values
+    out_y, out_x = base.shape[1], base.shape[2]
+    src = image.isel(x=slice(px0, px1), y=slice(py0, py1))
+    src = src.compute() if hasattr(src.data, "compute") else src
+    # Canvas spans the sliced source's own pixel extent, so datashader places every source pixel; the
+    # aggregated grid then lines up with ``base`` (same window, same resolution).
+    cvs = ds.Canvas(plot_width=out_x, plot_height=out_y, x_range=(px0, px1), y_range=(py0, py1))
+    agg = np.stack(
+        [np.asarray(cvs.raster(src.isel(c=i), downsample_method=downsample_method).values) for i in range(src.sizes["c"])],
+        axis=0,
+    )
+    base.values = agg.astype(base.dtype, copy=False)
+    return base
+
+
 def _rasterize_if_necessary_datashader(
     image: DataArray,
     dpi: float,
@@ -886,12 +949,18 @@ def _rasterize_if_necessary_datashader(
     coordinate_system: str,
     extent: dict[str, tuple[float, float]],
     downsample_method: str,
+    crop: BBox | None = None,
 ) -> DataArray:
     """Downsample to canvas resolution with a configurable datashader reduction.
 
     Used by ``render_images(method='datashader')`` so sparse images (mostly
     zeros, rare non-zero pixels) survive the downsample step instead of
     being averaged away by the default mean aggregation.
+
+    When ``crop`` (a cs-space ``(x0, y0, x1, y1)`` box) is given, the aggregation is restricted to the
+    window: ``rasterize`` supplies the correct window coords/transform and only the sliced source is
+    materialized (so the full image is never read — this scales to Visium HD). Falls back to the full
+    path for rotation/shear or an empty window.
     """
     has_c_dim = len(image.shape) == 3
     y_dims, x_dims = (image.shape[1], image.shape[2]) if has_c_dim else image.shape
@@ -899,6 +968,12 @@ def _rasterize_if_necessary_datashader(
     target_y_dims = int(dpi * height)
     target_x_dims = int(dpi * width)
 
+    if crop is not None:
+        windowed = _datashader_window_image(image, crop, coordinate_system, target_x_dims, target_y_dims, downsample_method)
+        if windowed is not None:
+            return windowed
+        # rotation/shear or empty window: fall through to the full render (axis limits clip)
+
     if y_dims <= target_y_dims and x_dims <= target_x_dims:
         return image
 
@@ -940,11 +1015,16 @@ def _multiscale_to_spatial_image(
     height: float,
     scale: str | None = None,
     is_label: bool = False,
+    crop: BBox | None = None,
+    extent: dict[str, tuple[float, float]] | None = None,
 ) -> DataArray:
     """Extract the DataArray to be rendered from a multiscale image.
 
     From the `DataTree`, the scale that fits the given image size and dpi most is selected
     and returned. In case the lowest resolution is still too high, a rasterization step is added.
+    When ``crop`` (a cs-space ``(x0, y0, x1, y1)`` box) and the full ``extent`` are given, only the
+    crop window is shown, so the level is chosen so the *window* reaches the target resolution rather
+    than the full image — otherwise a Visium-HD crop would upsample a coarse pyramid level.
 
     Parameters
     ----------
@@ -987,6 +1067,16 @@ def _multiscale_to_spatial_image(
         optimal_x = width * dpi
         optimal_y = height * dpi
 
+        if crop is not None and extent is not None:
+            # Only the crop window is drawn, so the full image needs enough pixels that the window
+            # portion reaches the target — scale the target up by full/window per axis.
+            full_x = float(extent["x"][1]) - float(extent["x"][0])
+            full_y = float(extent["y"][1]) - float(extent["y"][0])
+            win_x, win_y = crop.x1 - crop.x0, crop.y1 - crop.y0
+            if win_x > 0 and win_y > 0:
+                optimal_x *= full_x / win_x
+                optimal_y *= full_y / win_y
+
         # Pick the lowest-resolution scale where both x and y are >= the
         # target pixel count.  Falls back to highest available resolution.
         optimal_scale = scales[-1]
@@ -1459,6 +1549,62 @@ def _element_extent_fast(
     return {"x": (float(tc[:, 0].min()), float(tc[:, 0].max())), "y": (float(tc[:, 1].min()), float(tc[:, 1].max()))}
 
 
+# --- Fast bbox subset for on-the-fly cropping (points/shapes only) ---------------------------------
+# Drop rows/geometries outside a coordinate-system bbox before the expensive draw, without rebuilding
+# a SpatialData object or re-joining tables (which is what makes spatialdata.bounding_box_query slow).
+# Images/labels are not subset here: slicing a DataArray does not update what get_extent/imshow read,
+# so a sliced image is drawn in the wrong place; they rely on axis-limit clipping instead.
+
+
+def _bbox_to_element_space(element: Any, coordinate_system: str, bbox: BBox) -> BBox | None:
+    """Map a coordinate-system ``bbox`` to element-native coords.
+
+    Inverts the axis-aligned affine on the box corners and re-sorts so the result has ``x0 <= x1`` and
+    ``y0 <= y1`` even when the transform flips an axis (a negative scale would otherwise yield a
+    reversed box that silently empties ``.cx``/``.sel``). Returns ``None`` for rotation/shear, where a
+    box does not map to a box; the caller then skips the fast subset.
+    """
+    matrix = get_transformation(element, get_all=True)[coordinate_system].to_affine_matrix(("x", "y"), ("x", "y"))
+    affine = matrix[:2, :2]
+    if not _is_axis_aligned(affine):
+        return None
+    x0, y0, x1, y1 = bbox
+    corners = np.array([[x0, y0], [x1, y0], [x0, y1], [x1, y1]], dtype=float)
+    elem = (corners - matrix[:2, 2]) @ np.linalg.inv(affine).T
+    return BBox(float(elem[:, 0].min()), float(elem[:, 1].min()), float(elem[:, 0].max()), float(elem[:, 1].max()))
+
+
+def _bbox_mask_points(x: ArrayLike, y: ArrayLike, elem_bbox: BBox) -> ArrayLike:
+    """Boolean keep-mask for points whose element-native ``(x, y)`` fall inside ``elem_bbox``."""
+    x0, y0, x1, y1 = elem_bbox
+    x, y = np.asarray(x), np.asarray(y)
+    return (x >= x0) & (x <= x1) & (y >= y0) & (y <= y1)
+
+
+def _bbox_mask_shapes(shapes: GeoDataFrame, elem_bbox: BBox) -> ArrayLike:
+    """Boolean keep-mask for shapes whose bounds intersect ``elem_bbox`` (element-native coords).
+
+    Bbox-intersect, not clip: a boundary-crossing shape is kept whole (axis-limit clipping trims it
+    visually). Circles are stored as ``Point`` + ``radius``, so a centroid-only test would drop circles
+    whose body overlaps the box. Rows are handled per-geometry (a frame may mix circles and polygons):
+    polygons use the geometry's own bounds; circle rows expand their degenerate ``Point`` bounds by the
+    radius. This is the masking form of geopandas ``.cx`` (which returns a frame, not a mask we can align
+    to ``color_spec``).
+    """
+    x0, y0, x1, y1 = elem_bbox
+    geom = shapes.geometry
+    b = geom.bounds  # C-level; columns minx, miny, maxx, maxy (degenerate minx==maxx==cx for Points)
+    minx, miny, maxx, maxy = (b[c].to_numpy() for c in ("minx", "miny", "maxx", "maxy"))
+    is_point = (geom.geom_type == "Point").to_numpy()
+    if is_point.any() and ShapesModel.RADIUS_KEY in shapes:  # expand circle rows by their radius
+        r = np.asarray(shapes[ShapesModel.RADIUS_KEY], dtype=float)
+        minx = np.where(is_point, minx - r, minx)
+        miny = np.where(is_point, miny - r, miny)
+        maxx = np.where(is_point, maxx + r, maxx)
+        maxy = np.where(is_point, maxy + r, maxy)
+    return (maxx >= x0) & (minx <= x1) & (maxy >= y0) & (miny <= y1)
+
+
 def _fast_extent(element: Any, coordinate_system: str) -> dict[str, tuple[float, float]]:
     """Element extent via the fast corner-transform.
 
diff --git a/tests/_images/Show_crop_circles.png b/tests/_images/Show_crop_circles.png
new file mode 100644
index 00000000..fcd08b1a
Binary files /dev/null and b/tests/_images/Show_crop_circles.png differ
diff --git a/tests/_images/Show_crop_image.png b/tests/_images/Show_crop_image.png
new file mode 100644
index 00000000..59c50cbb
Binary files /dev/null and b/tests/_images/Show_crop_image.png differ
diff --git a/tests/_images/Show_crop_labels.png b/tests/_images/Show_crop_labels.png
new file mode 100644
index 00000000..6c935ee3
Binary files /dev/null and b/tests/_images/Show_crop_labels.png differ
diff --git a/tests/_images/Show_crop_layered_elements.png b/tests/_images/Show_crop_layered_elements.png
new file mode 100644
index 00000000..663dbbe2
Binary files /dev/null and b/tests/_images/Show_crop_layered_elements.png differ
diff --git a/tests/_images/Show_crop_points.png b/tests/_images/Show_crop_points.png
new file mode 100644
index 00000000..e88eb9eb
Binary files /dev/null and b/tests/_images/Show_crop_points.png differ
diff --git a/tests/_images/Show_crop_shapes.png b/tests/_images/Show_crop_shapes.png
new file mode 100644
index 00000000..655f2019
Binary files /dev/null and b/tests/_images/Show_crop_shapes.png differ
diff --git a/tests/pl/test_show.py b/tests/pl/test_show.py
index 1582560b..76dc3d8d 100644
--- a/tests/pl/test_show.py
+++ b/tests/pl/test_show.py
@@ -3,10 +3,13 @@
 
 import matplotlib
 import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
 import pytest
 import scanpy as sc
 from matplotlib.figure import Figure
 from spatialdata import SpatialData
+from spatialdata.models import PointsModel
 from spatialdata.transformations import Identity, set_transformation
 
 import spatialdata_plot  # noqa: F401
@@ -29,6 +32,34 @@ class TestShow(PlotTester, metaclass=PlotTesterMeta):
     def test_plot_pad_extent_adds_padding(self, sdata_blobs: SpatialData):
         sdata_blobs.pl.render_images(element="blobs_image").pl.show(pad_extent=100)
 
+    def test_plot_crop_image(self, sdata_blobs: SpatialData):
+        """Visual test: crop_coord windows an image to the box (#764)."""
+        sdata_blobs.pl.render_images("blobs_image").pl.show(crop_coord=(150, 400, 150, 400))
+
+    def test_plot_crop_points(self, sdata_blobs: SpatialData):
+        """Visual test: crop_coord subsets points to the box, colours from the full element (#764)."""
+        sdata_blobs.pl.render_points("blobs_points", color="genes", size=20).pl.show(crop_coord=(150, 400, 150, 400))
+
+    def test_plot_crop_shapes(self, sdata_blobs: SpatialData):
+        """Visual test: crop_coord subsets polygons to the box (#764)."""
+        sdata_blobs.pl.render_shapes("blobs_polygons").pl.show(crop_coord=(150, 400, 150, 400))
+
+    def test_plot_crop_circles(self, sdata_blobs: SpatialData):
+        """Visual test: crop_coord keeps circles whose body overlaps the box (radius-aware, #764)."""
+        sdata_blobs.pl.render_shapes("blobs_circles").pl.show(crop_coord=(150, 400, 150, 400))
+
+    def test_plot_crop_labels(self, sdata_blobs: SpatialData):
+        """Visual test: crop_coord windows a labels layer to the box (#764)."""
+        sdata_blobs.pl.render_labels("blobs_labels", color="channel_0_sum").pl.show(crop_coord=(150, 400, 150, 400))
+
+    def test_plot_crop_layered_elements(self, sdata_blobs: SpatialData):
+        """Visual test: layered image + labels both clip to the same crop box (#764)."""
+        (
+            sdata_blobs.pl.render_images("blobs_image")
+            .pl.render_labels("blobs_labels", fill_alpha=0.5)
+            .pl.show(crop_coord=(150, 400, 150, 400))
+        )
+
     def test_plot_xlabel_ylabel(self, sdata_blobs: SpatialData):
         """Visual test: xlabel/ylabel label the axes (feature for #763)."""
         sdata_blobs.pl.render_images(element="blobs_image").pl.show(xlabel="x (µm)", ylabel="y (µm)")
@@ -142,6 +173,264 @@ def test_title_empty_string_suppresses_title(self, sdata_blobs: SpatialData):
         plt.close("all")
 
 
+def test_crop_sets_exact_axis_limits(sdata_blobs: SpatialData):
+    """crop_coord=(xmin, xmax, ymin, ymax) pins the view to the box; y is inverted (top-left origin)."""
+    ax = sdata_blobs.pl.render_points().pl.show(crop_coord=(100, 300, 120, 260), return_ax=True, show=False)
+    assert ax.get_xlim() == pytest.approx((100, 300))
+    assert ax.get_ylim() == pytest.approx((260, 120))  # set_ylim(ymax, ymin)
+    plt.close("all")
+
+
+def test_crop_ignores_pad_extent(sdata_blobs: SpatialData):
+    """pad_extent must not widen a crop box (the view is exactly the box)."""
+    ax = sdata_blobs.pl.render_points().pl.show(crop_coord=(100, 300, 120, 260), pad_extent=50, return_ax=True, show=False)
+    assert ax.get_xlim() == pytest.approx((100, 300))
+    assert ax.get_ylim() == pytest.approx((260, 120))
+    plt.close("all")
+
+
+def test_crop_reduces_points_drawn(sdata_blobs: SpatialData):
+    """The fast subset draws fewer points than the full render."""
+
+    def n_offsets(ax):
+        return sum(len(c.get_offsets()) for c in ax.collections if hasattr(c, "get_offsets"))
+
+    full = sdata_blobs.pl.render_points().pl.show(return_ax=True, show=False)
+    cropped = sdata_blobs.pl.render_points().pl.show(crop_coord=(100, 300, 120, 260), return_ax=True, show=False)
+    assert 0 < n_offsets(cropped) < n_offsets(full)
+    plt.close("all")
+
+
+def test_crop_continuous_color_domain_from_full_element():
+    """Auto-scaled color range must come from the full element, so cropped colors match uncropped."""
+    rng = np.random.default_rng(0)
+    df = pd.DataFrame({"x": rng.uniform(0, 100, 2000), "y": rng.uniform(0, 100, 2000), "val": rng.uniform(0, 1, 2000)})
+    df.loc[0, ["x", "y", "val"]] = [90, 90, 10.0]  # an extreme value far outside the crop box
+    sdata = SpatialData(points={"p": PointsModel.parse(df, transformations={"global": Identity()})})
+
+    def vrange(ax):
+        for c in ax.collections:
+            if getattr(c, "norm", None) is not None and c.norm.vmax is not None:
+                return (c.norm.vmin, c.norm.vmax)
+        return None
+
+    full = sdata.pl.render_points("p", color="val").pl.show(return_ax=True, show=False)
+    cropped = sdata.pl.render_points("p", color="val").pl.show(crop_coord=(20, 50, 30, 60), return_ax=True, show=False)
+    assert vrange(cropped) == pytest.approx(vrange(full))
+    plt.close("all")
+
+
+def test_crop_transfunc_norm_matches_uncropped():
+    """The pinned norm must use the transfunc'd full-element range, so crop+transfunc matches uncropped."""
+    rng = np.random.default_rng(0)
+    df = pd.DataFrame({"x": rng.uniform(0, 100, 3000), "y": rng.uniform(0, 100, 3000), "val": rng.uniform(0, 1, 3000)})
+    df.loc[0, ["x", "y", "val"]] = [90, 90, 10.0]  # extreme value outside the crop box
+    sdata = SpatialData(points={"p": PointsModel.parse(df, transformations={"global": Identity()})})
+
+    def vrange(ax):
+        for c in ax.collections:
+            if getattr(c, "norm", None) is not None and c.norm.vmax is not None:
+                return (c.norm.vmin, c.norm.vmax)
+        return None
+
+    full = sdata.pl.render_points("p", color="val", transfunc=np.log1p, method="matplotlib").pl.show(
+        return_ax=True, show=False
+    )
+    cropped = sdata.pl.render_points("p", color="val", transfunc=np.log1p, method="matplotlib").pl.show(
+        crop_coord=(20, 50, 30, 60), return_ax=True, show=False
+    )
+    assert vrange(cropped) == pytest.approx(vrange(full))  # log1p range, not the raw range
+    plt.close("all")
+
+
+def test_crop_datashader_autoscales_over_window():
+    """Datashader crop autoscales over the visible window: a value far outside the box can't recolor it."""
+    rng = np.random.default_rng(0)
+    base = pd.DataFrame({"x": rng.uniform(20, 50, 12000), "y": rng.uniform(30, 60, 12000), "val": rng.uniform(0, 1, 12000)})
+    outside = pd.DataFrame({"x": [200.0], "y": [200.0], "val": [1000.0]})  # far outside the crop window
+    s_a = SpatialData(points={"p": PointsModel.parse(base, transformations={"global": Identity()})})
+    s_b = SpatialData(
+        points={"p": PointsModel.parse(pd.concat([base, outside], ignore_index=True), transformations={"global": Identity()})}
+    )
+
+    def raster(s):
+        ax = s.pl.render_points("p", color="val", method="datashader").pl.show(
+            crop_coord=(20, 50, 30, 60), return_ax=True, show=False
+        )
+        (im,) = ax.get_images()
+        arr = np.asarray(im.get_array()).copy()
+        plt.close("all")
+        return arr
+
+    # If the norm leaked the full-data range (old bug), the extreme value would squish the window's colors.
+    np.testing.assert_array_equal(raster(s_a), raster(s_b))
+
+
+def test_crop_multiscale_selects_finer_level():
+    """A crop must pick a pyramid level fine enough for the WINDOW, not the whole image (Visium HD)."""
+    from spatialdata.models import Image2DModel
+
+    from spatialdata_plot.pl.render_params import BBox
+    from spatialdata_plot.pl.utils import _multiscale_to_spatial_image
+
+    n = 800
+    rng = np.random.default_rng(0)
+    tree = Image2DModel.parse(rng.random((1, n, n), dtype=np.float32), dims=("c", "y", "x"), scale_factors=[2, 2])
+    extent = {"x": (0.0, float(n)), "y": (0.0, float(n))}
+    coarse = _multiscale_to_spatial_image(tree, dpi=10, width=5, height=5)  # target ~50px over the full image
+    fine = _multiscale_to_spatial_image(
+        tree, dpi=10, width=5, height=5, crop=BBox(0.0, 0.0, 80.0, 80.0), extent=extent  # 10% window -> 10x boost
+    )
+    assert fine.shape[-1] > coarse.shape[-1]
+
+
+def test_crop_image_rasterizes_only_window():
+    """A cropped large image is rasterized to the window at figure resolution, not the full image then clipped.
+
+    Regression for #764: rasterize() maps the crop bbox through the element transform, so placement is
+    correct even under a Scale+Translation (which a naive .sel would mis-place) and only the window is read.
+    """
+    from spatialdata.models import Image2DModel
+    from spatialdata.transformations import Scale, Sequence, Translation
+
+    n = 3000
+    rng = np.random.default_rng(0)
+    transform = Sequence([Scale([2.0, 2.0], axes=("x", "y")), Translation([1000.0, 500.0], axes=("x", "y"))])
+    img = Image2DModel.parse(
+        rng.random((1, n, n), dtype=np.float32), dims=("c", "y", "x"), transformations={"global": transform}
+    )
+    sdata = SpatialData(images={"img": img})
+    # full world extent x=(1000, 7000), y=(500, 6500); crop a small central window
+    ax = sdata.pl.render_images("img").pl.show(crop_coord=(3500, 3900, 3500, 3900), return_ax=True, show=False)
+
+    assert ax.get_xlim() == pytest.approx((3500, 3900))
+    assert ax.get_ylim() == pytest.approx((3900, 3500))  # inverted y
+    # the rendered raster covers the window at ~figure resolution, not the 3000-px source
+    (im,) = ax.get_images()
+    assert max(im.get_array().shape[:2]) < n // 2
+    plt.close("all")
+
+
+def test_crop_datashader_image_rasterizes_only_window():
+    """method='datashader' images are windowed under crop too, not full-rendered then clipped (#764 F4)."""
+    from spatialdata.models import Image2DModel
+    from spatialdata.transformations import Scale, Sequence, Translation
+
+    n = 3000
+    rng = np.random.default_rng(0)
+    transform = Sequence([Scale([2.0, 2.0], axes=("x", "y")), Translation([1000.0, 500.0], axes=("x", "y"))])
+    img = Image2DModel.parse(
+        rng.random((1, n, n), dtype=np.float32), dims=("c", "y", "x"), transformations={"global": transform}
+    )
+    sdata = SpatialData(images={"img": img})
+    ax = sdata.pl.render_images("img", method="datashader").pl.show(
+        crop_coord=(3500, 3900, 3500, 3900), return_ax=True, show=False
+    )
+    assert ax.get_xlim() == pytest.approx((3500, 3900))
+    assert ax.get_ylim() == pytest.approx((3900, 3500))  # inverted y
+    (im,) = ax.get_images()
+    assert max(im.get_array().shape[:2]) < n // 2  # window at figure resolution, not the full source
+    plt.close("all")
+
+
+def _grid_labels_sdata():
+    """A 2000x2000 label raster of 400 block-instances with a Scale+Translation, plus a table with a
+    categorical and a plain-string colour column. Large enough that rasterize()/windowing engages."""
+    import anndata as ad
+    from spatialdata.models import Labels2DModel, TableModel
+    from spatialdata.transformations import Scale, Sequence, Translation
+
+    n = 2000
+    rng = np.random.default_rng(1)
+    lab = np.zeros((n, n), dtype=np.int32)
+    k = 0
+    for i in range(20):
+        for j in range(20):
+            k += 1
+            lab[i * 100 : (i + 1) * 100, j * 100 : (j + 1) * 100] = k
+    transform = Sequence([Scale([2.0, 2.0], axes=("x", "y")), Translation([1000.0, 500.0], axes=("x", "y"))])
+    labels = Labels2DModel.parse(lab, dims=("y", "x"), transformations={"global": transform})
+    obs = pd.DataFrame(
+        {
+            "instance_id": np.arange(1, 401),
+            "region": pd.Categorical(["labels"] * 400),
+            "ct_cat": pd.Categorical(rng.choice(list("ABCDE"), size=400), categories=list("ABCDE")),
+            "ct_str": rng.choice(list("ABCDE"), size=400).astype(object),
+        }
+    )
+    table = TableModel.parse(ad.AnnData(obs=obs), region="labels", region_key="region", instance_key="instance_id")
+    return SpatialData(labels={"labels": labels}, tables={"table": table})
+
+
+def _legend_colors(ax):
+    leg = ax.get_legend()
+    out = {}
+    if leg is None:
+        return out
+    for text, handle in zip(leg.get_texts(), leg.legend_handles):
+        for attr in ("get_facecolor", "get_color"):
+            try:
+                v = np.ravel(getattr(handle, attr)())
+                if v.size >= 3:
+                    out[text.get_text()] = tuple(np.round(v[:3], 3))
+                    break
+            except (AttributeError, TypeError):
+                pass
+    return out
+
+
+def test_crop_labels_placement_and_empty_window():
+    """Labels crop pins the view to the box; a window with no labels renders without crashing."""
+    sdata = _grid_labels_sdata()
+    ax = sdata.pl.render_labels("labels", color="ct_cat").pl.show(
+        crop_coord=(2600, 3000, 2200, 2600), return_ax=True, show=False
+    )
+    assert ax.get_xlim() == pytest.approx((2600, 3000))
+    assert ax.get_ylim() == pytest.approx((2600, 2200))  # inverted y
+    plt.close("all")
+    # a box far outside the data must not raise (empty-window guard)
+    sdata.pl.render_labels("labels", color="ct_cat").pl.show(crop_coord=(99000, 99400, 99000, 99400), show=False)
+    plt.close("all")
+
+
+@pytest.mark.parametrize("col", ["ct_cat", "ct_str"])
+def test_crop_labels_no_color_reshuffle(col):
+    """Cropped label colours must match the uncropped plot for shared categories (windowed dtype-Categorical;
+    plain-string falls back to full render so it stays stable too)."""
+    sdata = _grid_labels_sdata()
+    full = sdata.pl.render_labels("labels", color=col).pl.show(return_ax=True, show=False)
+    full_colors = _legend_colors(full)
+    plt.close("all")
+    cropped = sdata.pl.render_labels("labels", color=col).pl.show(
+        crop_coord=(2600, 3000, 2200, 2600), return_ax=True, show=False
+    )
+    crop_colors = _legend_colors(cropped)
+    plt.close("all")
+    shared = [c for c in set(full_colors) & set(crop_colors)]
+    assert shared  # the window keeps several categories
+    for c in shared:
+        assert full_colors[c] == pytest.approx(crop_colors[c], abs=0.02)
+
+
+def test_crop_invalid_order_raises(sdata_blobs: SpatialData):
+    with pytest.raises(ValueError, match="xmin < xmax and ymin < ymax"):
+        sdata_blobs.pl.render_points().pl.show(crop_coord=(300, 100, 120, 260), show=False)
+
+
+def test_crop_wrong_length_raises(sdata_blobs: SpatialData):
+    with pytest.raises(TypeError, match="tuple of four numbers"):
+        sdata_blobs.pl.render_points().pl.show(crop_coord=(100, 300, 120), show=False)
+
+
+def test_crop_multiple_coordinate_systems_raises(sdata_blobs: SpatialData):
+    """crop is one box in one CS's units; rendering several CS at once is rejected."""
+    set_transformation(sdata_blobs["blobs_points"], Identity(), to_coordinate_system="other")
+    with pytest.raises(ValueError, match="single coordinate system"):
+        sdata_blobs.pl.render_points().pl.show(
+            coordinate_systems=["global", "other"], crop_coord=(100, 300, 120, 260), show=False
+        )
+
+
 def test_fig_parameter_emits_deprecation_warning(sdata_blobs: SpatialData):
     """Passing fig= should emit a DeprecationWarning (regression for #204)."""
     fig = Figure()
diff --git a/tests/pl/test_utils.py b/tests/pl/test_utils.py
index fac573a5..7eeb4bb4 100644
--- a/tests/pl/test_utils.py
+++ b/tests/pl/test_utils.py
@@ -7,9 +7,10 @@
 import scanpy as sc
 import xarray as xr
 from anndata import AnnData
-from shapely.geometry import Point
+from shapely.geometry import Point, Polygon
 from spatialdata import SpatialData, get_centroids
 from spatialdata.models import Labels2DModel, PointsModel, ShapesModel, TableModel
+from spatialdata.transformations import Affine, Identity, Scale, Sequence, Translation
 
 import spatialdata_plot
 from spatialdata_plot.pl import measure_obs
@@ -19,7 +20,12 @@
     _datashader_map_aggregate_to_color,
 )
 from spatialdata_plot.pl.render_params import CmapParams, Color, ColorLike, colormap_with_alpha
-from spatialdata_plot.pl.utils import set_zero_in_cmap_to_transparent
+from spatialdata_plot.pl.utils import (
+    _bbox_mask_points,
+    _bbox_mask_shapes,
+    _bbox_to_element_space,
+    set_zero_in_cmap_to_transparent,
+)
 from tests.conftest import DPI, PlotTester, PlotTesterMeta
 
 sc.pl.set_rcParams_defaults()
@@ -1355,3 +1361,124 @@ def test_datashader_points_image_aligns_with_points_extent():
     assert (min(x0, x1), max(x0, x1)) == pytest.approx((0.1, 0.9))
     assert (min(y0, y1), max(y0, y1)) == pytest.approx((0.1, 0.9))
     plt.close("all")
+
+
+# --- bbox-subset helpers (on-the-fly crop) --------------------------------------------------------
+
+
+def _points(x, y):
+    return PointsModel.parse(pd.DataFrame({"x": x, "y": y}), transformations={"g": Identity()})
+
+
+def test_bbox_to_element_space_identity_round_trips():
+    pts = _points([0.0, 1.0], [0.0, 1.0])
+    assert _bbox_to_element_space(pts, "g", (2.0, 3.0, 8.0, 9.0)) == pytest.approx((2.0, 3.0, 8.0, 9.0))
+
+
+def test_bbox_to_element_space_inverts_scale_and_translation():
+    t = Sequence([Scale([2, 2], axes=("x", "y")), Translation([100, 50], axes=("x", "y"))])
+    pts = PointsModel.parse(pd.DataFrame({"x": [0.0], "y": [0.0]}), transformations={"g": t})
+    # cs box x[110,126] y[110,116] -> element x[5,13] y[30,33]
+    assert _bbox_to_element_space(pts, "g", (110.0, 110.0, 126.0, 116.0)) == pytest.approx((5.0, 30.0, 13.0, 33.0))
+
+
+def test_bbox_to_element_space_resorts_under_axis_flip():
+    # a negative x-scale flips the axis; the element box must keep x0 < x1 so .sel/mask don't empty
+    flip = Scale([-1, 1], axes=("x", "y"))
+    pts = PointsModel.parse(pd.DataFrame({"x": [0.0], "y": [0.0]}), transformations={"g": flip})
+    x0, y0, x1, y1 = _bbox_to_element_space(pts, "g", (2.0, 1.0, 6.0, 3.0))
+    assert x0 < x1 and y0 < y1
+
+
+def test_bbox_to_element_space_none_for_rotation():
+    theta = 0.3
+    rot = Affine(
+        [[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]],
+        input_axes=("x", "y"),
+        output_axes=("x", "y"),
+    )
+    pts = PointsModel.parse(pd.DataFrame({"x": [0.0], "y": [0.0]}), transformations={"g": rot})
+    assert _bbox_to_element_space(pts, "g", (0.0, 0.0, 1.0, 1.0)) is None
+
+
+def test_bbox_mask_points_matches_numpy_predicate():
+    x = np.array([1.0, 5.0, 9.0, 15.0])
+    y = np.array([1.0, 5.0, 9.0, 15.0])
+    mask = _bbox_mask_points(x, y, (2.0, 2.0, 8.0, 8.0))
+    np.testing.assert_array_equal(mask, [False, True, False, False])
+
+
+def test_bbox_mask_shapes_keeps_circle_near_edge():
+    # circle centred at (0,0) r=3: centroid outside box (1..2) but body overlaps -> must be kept
+    circ = ShapesModel.parse(
+        gpd.GeoDataFrame({"geometry": [Point(0, 0)], "radius": [3.0]}), transformations={"g": Identity()}
+    )
+    assert _bbox_mask_shapes(circ, (1.0, 1.0, 2.0, 2.0)).tolist() == [True]
+    # a centroid-only test (geopandas .cx) would have dropped it
+    assert len(circ.cx[1:2, 1:2]) == 0
+
+
+def test_bbox_mask_shapes_mixed_circle_and_polygon():
+    # a frame mixing a circle (Point+radius) and a polygon must stay radius-aware per row:
+    # the r=3 circle body reaches the box; the distant polygon does not.
+    gdf = gpd.GeoDataFrame(
+        {"geometry": [Point(0, 0), Polygon([(100, 100), (101, 100), (101, 101)])], "radius": [3.0, np.nan]}
+    )
+    shapes = ShapesModel.parse(gdf, transformations={"g": Identity()})
+    np.testing.assert_array_equal(_bbox_mask_shapes(shapes, (1.0, 1.0, 2.0, 2.0)), [True, False])
+
+
+def test_bbox_mask_shapes_polygon_intersect_and_empty():
+    polys = [Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]), Polygon([(100, 100), (110, 100), (110, 110), (100, 110)])]
+    gdf = ShapesModel.parse(gpd.GeoDataFrame({"geometry": polys}), transformations={"g": Identity()})
+    np.testing.assert_array_equal(_bbox_mask_shapes(gdf, (1.0, 1.0, 5.0, 5.0)), [True, False])
+    np.testing.assert_array_equal(_bbox_mask_shapes(gdf, (200.0, 200.0, 210.0, 210.0)), [False, False])
+
+
+# --- datashader image windowing (crop) ------------------------------------------------------------
+
+
+def _bright_block_image(n=200):
+    """A single-channel n×n image (Identity transform) with a bright block at y∈[120,140), x∈[40,60)."""
+    from spatialdata.models import Image2DModel
+
+    arr = np.zeros((1, n, n), dtype=np.float32)
+    arr[0, 120:140, 40:60] = 1.0
+    return Image2DModel.parse(arr, dims=("c", "y", "x"), transformations={"g": Identity()})
+
+
+def test_datashader_window_image_places_block_at_correct_relative_position():
+    """The windowed aggregate must land the bright block where the crop box puts it (alignment fix)."""
+    from spatialdata_plot.pl.render_params import BBox
+    from spatialdata_plot.pl.utils import _datashader_window_image
+
+    img = _bright_block_image()
+    # window x[20,80] y[100,160]; block centre world (x=50, y=130) -> rel (0.5, 0.5) -> output centre
+    out = _datashader_window_image(img, BBox(20, 100, 80, 160), "g", 60, 60, "max")
+    assert out is not None
+    a = np.asarray(out.values[0])
+    ys, xs = np.where(a > 0.5)
+    assert ys.size > 0
+    oy, ox = a.shape
+    assert 0.3 * oy <= ys.mean() <= 0.7 * oy  # near vertical centre, not shifted to an edge
+    assert 0.3 * ox <= xs.mean() <= 0.7 * ox
+
+
+def test_datashader_window_image_clamps_at_image_edge():
+    """A window straddling the image edge still returns a raster (source is clamped, not dropped)."""
+    from spatialdata_plot.pl.render_params import BBox
+    from spatialdata_plot.pl.utils import _datashader_window_image
+
+    img = _bright_block_image()
+    out = _datashader_window_image(img, BBox(-50, 100, 50, 160), "g", 60, 60, "max")  # left half outside
+    assert out is not None
+    assert np.asarray(out.values[0]).max() > 0.5  # the in-bounds part of the block survives
+
+
+def test_datashader_window_image_none_when_fully_outside():
+    """A window with no source pixels falls back (returns None) instead of building an empty canvas."""
+    from spatialdata_plot.pl.render_params import BBox
+    from spatialdata_plot.pl.utils import _datashader_window_image
+
+    img = _bright_block_image()
+    assert _datashader_window_image(img, BBox(500, 500, 600, 600), "g", 60, 60, "max") is None