diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6e244f2e..f2ce9b59 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,20 +10,24 @@ and this project adheres to [Semantic Versioning][].
 
 ## [0.2.3] - tbd
 
+### Added
+
+-   Datashader support for points and shapes (#244)
+
 ### Changed
 
--   All parameters are now provided for a single element. If element in pl.render is None then this value will be broadcasted
+-   All parameters are now provided for a single element (#272)
 
 ### Fixed
 
 -   Fix color assignment for NaN values (#257)
--   Fix channel str support #221
+-   Zorder of rendering now strictly follows the order of the render_x calls (#244)
 
 ## [0.2.2] - 2024-05-02
 
 ### Fixed
 
--   Fixed `fill_alpha` ignoring `alpha` channel from custom cmap
+-   Fixed `fill_alpha` ignoring `alpha` channel from custom cmap (#236)
 -   Fix channel str support (#221)
 
 ## [0.2.1] - 2024-03-26
diff --git a/src/spatialdata_plot/pl/basic.py b/src/spatialdata_plot/pl/basic.py
index 5538c2a0..41025541 100644
--- a/src/spatialdata_plot/pl/basic.py
+++ b/src/spatialdata_plot/pl/basic.py
@@ -166,6 +166,7 @@ def render_shapes(
         cmap: Colormap | str | None = None,
         norm: bool | Normalize = False,
         scale: float | int = 1.0,
+        method: str | None = None,
         table_name: str | None = None,
         **kwargs: Any,
     ) -> sd.SpatialData:
@@ -219,6 +220,9 @@ def render_shapes(
             Colormap normalization for continuous annotations.
         scale : float | int, default 1.0
             Value to scale circles, if present.
+        method : str | None, optional
+            Whether to use 'matplotlib' and 'datashader'. When None, the method is
+            chosen based on the size of the data.
         table_name: str | None
             Name of the table containing the color(s) columns. If one name is given than the table is used for each
             spatial element to be plotted if the table annotates it. If you want to use different tables for particular
@@ -253,6 +257,7 @@ def render_shapes(
             norm=norm,
             scale=scale,
             table_name=table_name,
+            method=method,
         )
 
         sdata = self._copy()
@@ -281,6 +286,8 @@ def render_shapes(
                 fill_alpha=param_values["fill_alpha"],
                 transfunc=kwargs.get("transfunc", None),
                 table_name=param_values["table_name"],
+                zorder=n_steps,
+                method=method,
             )
             n_steps += 1
 
@@ -298,6 +305,7 @@ def render_points(
         cmap: Colormap | str | None = None,
         norm: None | Normalize = None,
         size: float | int = 1.0,
+        method: str | None = None,
         table_name: str | None = None,
         **kwargs: Any,
     ) -> sd.SpatialData:
@@ -342,6 +350,9 @@ def render_points(
             Colormap normalization for continuous annotations.
         size : float | int, default 1.0
             Size of the points
+        method : str | None, optional
+            Whether to use 'matplotlib' and 'datashader'. When None, the method is
+            chosen based on the size of the data.
         table_name: str | None
             Name of the table containing the color(s) columns. If one name is given than the table is used for each
             spatial element to be plotted if the table annotates it. If you want to use different tables for particular
@@ -368,6 +379,12 @@ def render_points(
             table_name=table_name,
         )
 
+        if method is not None:
+            if not isinstance(method, str):
+                raise TypeError("Parameter 'method' must be a string.")
+            if method not in ["matplotlib", "datashader"]:
+                raise ValueError("Parameter 'method' must be either 'matplotlib' or 'datashader'.")
+
         sdata = self._copy()
         sdata = _verify_plotting_tree(sdata)
         n_steps = len(sdata.plotting_tree.keys())
@@ -391,6 +408,8 @@ def render_points(
                 transfunc=kwargs.get("transfunc", None),
                 size=param_values["size"],
                 table_name=param_values["table_name"],
+                zorder=n_steps,
+                method=method,
             )
             n_steps += 1
 
@@ -504,6 +523,7 @@ def render_images(
                 alpha=param_values["alpha"],
                 percentiles_for_norm=param_values["percentiles_for_norm"],
                 scale=param_values["scale"],
+                zorder=n_steps,
             )
             n_steps += 1
 
@@ -515,14 +535,14 @@ def render_labels(
         element: str | None = None,
         color: str | None = None,
         groups: list[str] | str | None = None,
-        contour_px: int = 3,
+        contour_px: int | None = 3,
         outline: bool = False,
         palette: list[str] | str | None = None,
         cmap: Colormap | str | None = None,
         norm: Normalize | None = None,
         na_color: ColorLike | None = (0.0, 0.0, 0.0, 0.0),
         outline_alpha: float | int = 1.0,
-        fill_alpha: float | int = 0.3,
+        fill_alpha: float | int = 0.4,
         scale: str | None = None,
         table_name: str | None = None,
         **kwargs: Any,
@@ -626,6 +646,7 @@ def render_labels(
                 transfunc=kwargs.get("transfunc", None),
                 scale=param_values["scale"],
                 table_name=param_values["table_name"],
+                zorder=n_steps,
             )
             n_steps += 1
         return sdata
diff --git a/src/spatialdata_plot/pl/render.py b/src/spatialdata_plot/pl/render.py
index 09c85629..25b1dbf6 100644
--- a/src/spatialdata_plot/pl/render.py
+++ b/src/spatialdata_plot/pl/render.py
@@ -6,6 +6,7 @@
 from typing import Union
 
 import dask
+import datashader as ds
 import geopandas as gpd
 import matplotlib
 import matplotlib.transforms as mtransforms
@@ -15,12 +16,14 @@
 import spatialdata as sd
 from anndata import AnnData
 from datatree import DataTree
+from matplotlib.cm import ScalarMappable
 from matplotlib.colors import ListedColormap, Normalize
 from scanpy._settings import settings as sc_settings
 from spatialdata import get_extent
 from spatialdata.models import PointsModel, get_table_keys
 from spatialdata.transformations import (
     get_transformation,
+    set_transformation,
 )
 
 from spatialdata_plot._logging import logger
@@ -42,6 +45,7 @@
     _is_coercable_to_float,
     _map_color_seg,
     _maybe_set_colors,
+    _mpl_ax_contains_elements,
     _multiscale_to_spatial_image,
     _normalize,
     _rasterize_if_necessary,
@@ -70,7 +74,6 @@ def _render_shapes(
         filter_tables=bool(render_params.table_name),
     )
 
-    # for index, e in enumerate(elements):
     shapes = sdata[element]
 
     if (table_name := render_params.table_name) is None:
@@ -119,7 +122,6 @@ def _render_shapes(
         color_vector = [render_params.cmap_params.na_color]
 
     # filter by `groups`
-
     if isinstance(groups, list) and color_source_vector is not None:
         mask = color_source_vector.isin(groups)
         shapes = shapes[mask]
@@ -128,33 +130,6 @@ def _render_shapes(
         color_vector = color_vector[mask]
     shapes = gpd.GeoDataFrame(shapes, geometry="geometry")
 
-    _cax = _get_collection_shape(
-        shapes=shapes,
-        s=render_params.scale,
-        c=color_vector,
-        render_params=render_params,
-        rasterized=sc_settings._vector_friendly,
-        cmap=render_params.cmap_params.cmap,
-        norm=norm,
-        fill_alpha=render_params.fill_alpha,
-        outline_alpha=render_params.outline_alpha,
-        # **kwargs,
-    )
-
-    # Sets the limits of the colorbar to the values instead of [0, 1]
-    if not norm and not values_are_categorical:
-        _cax.set_clim(min(color_vector), max(color_vector))
-
-    cax = ax.add_collection(_cax)
-
-    # Apply the transformation to the PatchCollection's paths
-    trans = get_transformation(sdata_filt[element], get_all=True)[coordinate_system]
-    affine_trans = trans.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))
-    trans = mtransforms.Affine2D(matrix=affine_trans)
-
-    for path in _cax.get_paths():
-        path.vertices = trans.transform(path.vertices)
-
     # Using dict.fromkeys here since set returns in arbitrary order
     # remove the color of NaN values, else it might be assigned to a category
     # order of color in the palette should agree to order of occurence
@@ -171,6 +146,126 @@ def _render_shapes(
         # False if user specified color-like with 'color' parameter
         colorbar = False if col_for_color is None else legend_params.colorbar
 
+    # Apply the transformation to the PatchCollection's paths
+    trans = get_transformation(sdata_filt.shapes[element], get_all=True)[coordinate_system]
+    affine_trans = trans.to_affine_matrix(input_axes=("x", "y"), output_axes=("x", "y"))
+    trans = mtransforms.Affine2D(matrix=affine_trans)
+
+    shapes = gpd.GeoDataFrame(shapes, geometry="geometry")
+
+    # Determine which method to use for rendering
+    method = render_params.method
+    if method is None:
+        method = "datashader" if len(shapes) > 10000 else "matplotlib"
+    elif method not in ["matplotlib", "datashader"]:
+        raise ValueError("Method must be either 'matplotlib' or 'datashader'.")
+    logger.info(f"Using {method}")
+
+    if method == "datashader":
+        trans = mtransforms.Affine2D(matrix=affine_trans) + ax.transData
+
+        extent = get_extent(sdata.shapes[element])
+        x_ext = extent["x"][1]
+        y_ext = extent["y"][1]
+        x_range = [0, x_ext]
+        y_range = [0, y_ext]
+        # round because we need integers
+        plot_width = int(np.round(x_range[1] - x_range[0]))
+        plot_height = int(np.round(y_range[1] - y_range[0]))
+
+        cvs = ds.Canvas(plot_width=plot_width, plot_height=plot_height, x_range=x_range, y_range=y_range)
+
+        _geometry = shapes["geometry"]
+        is_point = _geometry.type == "Point"
+
+        # Handle circles encoded as points with radius
+        if is_point.any():
+            scale = shapes[is_point]["radius"] * render_params.scale
+            sdata_filt.shapes[element].loc[is_point, "geometry"] = _geometry[is_point].buffer(scale)
+
+        # in case we are coloring by a column in table
+        if col_for_color is not None and col_for_color not in sdata_filt.shapes[element].columns:
+            sdata_filt.shapes[element][col_for_color] = (
+                color_vector if color_source_vector is None else color_source_vector
+            )
+        # Render shapes with datashader
+        color_by_categorical = col_for_color is not None and color_source_vector is not None
+        aggregate_with_sum = None
+        if col_for_color is not None and (render_params.groups is None or len(render_params.groups) > 1):
+            if color_by_categorical:
+                agg = cvs.polygons(
+                    sdata_filt.shapes[element], geometry="geometry", agg=ds.by(col_for_color, ds.count())
+                )
+            else:
+                agg = cvs.polygons(sdata_filt.shapes[element], geometry="geometry", agg=ds.sum(column=col_for_color))
+                # save min and max values for drawing the colorbar
+                aggregate_with_sum = (agg.min(), agg.max())
+        else:
+            agg = cvs.polygons(sdata_filt.shapes[element], geometry="geometry", agg=ds.count())
+
+        color_key = (
+            [x[:-2] for x in color_vector.categories.values]
+            if (type(color_vector) is pd.core.arrays.categorical.Categorical)
+            and (len(color_vector.categories.values) > 1)
+            else None
+        )
+
+        ds_result = (
+            ds.tf.shade(
+                agg,
+                cmap=color_vector[0][:-2],
+                color_key=color_key,
+                min_alpha=np.min([150, render_params.fill_alpha * 255]),
+            )
+            if color_by_categorical or col_for_color is None
+            else ds.tf.shade(
+                agg,
+                cmap=render_params.cmap_params.cmap,
+            )
+        )
+        rgba_image = np.transpose(ds_result.to_numpy().base, (0, 1, 2))
+        _cax = ax.imshow(rgba_image, cmap=palette, zorder=render_params.zorder)
+        _cax.set_transform(trans)
+        cax = ax.add_image(_cax)
+        if aggregate_with_sum is not None:
+            cax = ScalarMappable(
+                norm=matplotlib.colors.Normalize(vmin=aggregate_with_sum[0], vmax=aggregate_with_sum[1]),
+                cmap=render_params.cmap_params.cmap,
+            )
+
+    elif method == "matplotlib":
+        _cax = _get_collection_shape(
+            shapes=shapes,
+            s=render_params.scale,
+            c=color_vector,
+            render_params=render_params,
+            rasterized=sc_settings._vector_friendly,
+            cmap=render_params.cmap_params.cmap,
+            norm=norm,
+            fill_alpha=render_params.fill_alpha,
+            outline_alpha=render_params.outline_alpha,
+            zorder=render_params.zorder,
+            # **kwargs,
+        )
+        cax = ax.add_collection(_cax)
+
+        # Transform the paths in PatchCollection
+        for path in _cax.get_paths():
+            path.vertices = trans.transform(path.vertices)
+            cax = ax.add_collection(_cax)
+
+    # Sets the limits of the colorbar to the values instead of [0, 1]
+    if not norm and not values_are_categorical:
+        _cax.set_clim(min(color_vector), max(color_vector))
+
+    if not (len(set(color_vector)) == 1 and list(set(color_vector))[0] == to_hex(render_params.cmap_params.na_color)):
+        # necessary in case different shapes elements are annotated with one table
+        if color_source_vector is not None and render_params.col_for_color is not None:
+            color_source_vector = color_source_vector.remove_unused_categories()
+
+        # False if user specified color-like with 'color' parameter
+        colorbar = False if render_params.col_for_color is None else legend_params.colorbar
+
         _ = _decorate_axs(
             ax=ax,
             cax=cax,
@@ -236,6 +331,8 @@ def _render_points(
 
     if groups is not None and col_for_color is not None:
         points = points[points[col_for_color].isin(groups)]
+        if len(points) <= 0:
+            raise ValueError(f"None of the groups {groups} could be found in the column '{col_for_color}'.")
 
     # we construct an anndata to hack the plotting functions
     if table_name is None:
@@ -251,11 +348,16 @@ def _render_points(
         )
         sdata_filt[table_name] = adata
 
-    # we can do this because of dealing with a copy
+    # we can modify the sdata because of dealing with a copy
 
     # Convert back to dask dataframe to modify sdata
+    transformation_in_cs = sdata_filt.points[element].attrs["transform"][coordinate_system]
     points = dask.dataframe.from_pandas(points, npartitions=1)
     sdata_filt.points[element] = PointsModel.parse(points, coordinates={"x": "x", "y": "y"})
+    # restore transformation in coordinate system of interest
+    set_transformation(
+        element=sdata_filt.points[element], transformation=transformation_in_cs, to_coordinate_system=coordinate_system
+    )
 
     if col_for_color is not None:
         cols = sc.get.obs_df(adata, col_for_color)
@@ -292,18 +394,98 @@ def _render_points(
     trans = mtransforms.Affine2D(matrix=affine_trans) + ax.transData
 
     norm = copy(render_params.cmap_params.norm)
-    _cax = ax.scatter(
-        adata[:, 0].X.flatten(),
-        adata[:, 1].X.flatten(),
-        s=render_params.size,
-        c=color_vector,
-        rasterized=sc_settings._vector_friendly,
-        cmap=render_params.cmap_params.cmap,
-        norm=norm,
-        alpha=render_params.alpha,
-        transform=trans,
-    )
-    cax = ax.add_collection(_cax)
+
+    method = render_params.method
+    if method is None:
+        method = "datashader" if len(points) > 10000 else "matplotlib"
+    elif method not in ["matplotlib", "datashader"]:
+        raise ValueError("Method must be either 'matplotlib' or 'datashader'.")
+
+    if method == "datashader":
+        # NOTE: s in matplotlib is in units of points**2
+        px = int(np.round(np.sqrt(render_params.size)))
+
+        extent = get_extent(sdata_filt.points[element], coordinate_system=coordinate_system)
+        x_ext = [min(0, extent["x"][0]), extent["x"][1]]
+        y_ext = [min(0, extent["y"][0]), extent["y"][1]]
+        previous_xlim = ax.get_xlim()
+        previous_ylim = ax.get_ylim()
+        # increase range if sth larger was rendered before
+        if _mpl_ax_contains_elements(ax):
+            x_ext = [min(x_ext[0], previous_xlim[0]), max(x_ext[1], previous_xlim[1])]
+            if ax.yaxis_inverted():  # case for e.g. images
+                y_ext = [min(y_ext[0], previous_ylim[1]), max(y_ext[1], previous_ylim[0])]
+            else:  # case for e.g. labels
+                y_ext = [min(y_ext[0], previous_ylim[0]), max(y_ext[1], previous_ylim[1])]
+        # round because we need integers
+        plot_width = int(np.round(x_ext[1] - x_ext[0]))
+        plot_height = int(np.round(y_ext[1] - y_ext[0]))
+
+        # use datashader for the visualization of points
+        cvs = ds.Canvas(plot_width=plot_width, plot_height=plot_height, x_range=x_ext, y_range=y_ext)
+
+        color_by_categorical = col_for_color is not None and points[col_for_color].values.dtype == object
+        aggregate_with_sum = None
+        if col_for_color is not None and (render_params.groups is None or len(render_params.groups) > 1):
+            if color_by_categorical:
+                agg = cvs.points(sdata_filt.points[element], "x", "y", agg=ds.by(col_for_color, ds.count()))
+            else:
+                agg = cvs.points(sdata_filt.points[element], "x", "y", agg=ds.sum(column=col_for_color))
+                # save min and max values for drawing the colorbar
+                aggregate_with_sum = (agg.min(), agg.max())
+        else:
+            agg = cvs.points(sdata_filt.points[element], "x", "y", agg=ds.count())
+
+        color_key = (
+            [x[:-2] for x in color_vector.categories.values]
+            if (type(color_vector) is pd.core.arrays.categorical.Categorical)
+            and (len(color_vector.categories.values) > 1)
+            else None
+        )
+        if color_by_categorical or col_for_color is None:
+            ds_result = ds.tf.shade(
+                ds.tf.spread(agg, px=px),
+                rescale_discrete_levels=True,
+                cmap=color_vector[0][:-2],
+                color_key=color_key,
+                min_alpha=np.min([150, render_params.alpha * 255]),  # value 150 is arbitrarily chosen
+            )
+        else:
+            ds_result = ds.tf.shade(
+                ds.tf.spread(agg, px=px),
+                rescale_discrete_levels=True,
+                cmap=render_params.cmap_params.cmap,
+            )
+
+        rbga_image = np.transpose(ds_result.to_numpy().base, (0, 1, 2))
+        cax = ax.imshow(rbga_image, zorder=render_params.zorder, alpha=render_params.alpha)
+        if aggregate_with_sum is not None:
+            cax = ScalarMappable(
+                norm=matplotlib.colors.Normalize(vmin=aggregate_with_sum[0], vmax=aggregate_with_sum[1]),
+                cmap=render_params.cmap_params.cmap,
+            )
+
+    elif method == "matplotlib":
+        # update axis limits if plot was empty before (necessary if datashader comes after)
+        update_parameters = not _mpl_ax_contains_elements(ax)
+        _cax = ax.scatter(
+            adata[:, 0].X.flatten(),
+            adata[:, 1].X.flatten(),
+            s=render_params.size,
+            c=color_vector,
+            rasterized=sc_settings._vector_friendly,
+            cmap=render_params.cmap_params.cmap,
+            norm=norm,
+            alpha=render_params.alpha,
+            transform=trans,
+            zorder=render_params.zorder,
+        )
+        cax = ax.add_collection(_cax)
+        if update_parameters:
+            # necessary if points are plotted with mpl first and then with datashader
+            extent = get_extent(sdata_filt.points[element], coordinate_system=coordinate_system)
+            ax.set_xbound(extent["x"])
+            ax.set_ybound(extent["y"])
 
     if len(set(color_vector)) != 1 or list(set(color_vector))[0] != to_hex(render_params.cmap_params.na_color):
         if color_source_vector is None:
@@ -414,14 +596,13 @@ def _render_images(
             _get_linear_colormap(palette, "k")[0]
             if isinstance(palette, list) and all(isinstance(p, str) for p in palette)
             else render_params.cmap_params.cmap
-            # render_params.cmap_params.cmap if render_params.palette is None else _get_linear_colormap(palette, "k")[0]
         )
 
         # Overwrite alpha in cmap: https://stackoverflow.com/a/10127675
         cmap._init()
         cmap._lut[:, -1] = render_params.alpha
 
-        _ax_show_and_transform(layer, trans_data, ax, cmap=cmap)
+        _ax_show_and_transform(layer, trans_data, ax, cmap=cmap, zorder=render_params.zorder)
 
     # 2) Image has any number of channels but 1
     else:
@@ -437,10 +618,12 @@ def _render_images(
                     clip=True,
                 )
 
-            if not isinstance(render_params.cmap_params, list) and render_params.cmap_params.norm:
-                layers[c] = render_params.cmap_params.norm(layers[c])
-            elif isinstance(render_params.cmap_params, list) and render_params.cmap_params[ch_index].norm:
-                layers[c] = render_params.cmap_params[ch_index].norm(layers[c])
+            if not isinstance(render_params.cmap_params, list):
+                if render_params.cmap_params.norm is not None:
+                    layers[c] = render_params.cmap_params.norm(layers[c])
+            else:
+                if render_params.cmap_params[ch_index].norm is not None:
+                    layers[c] = render_params.cmap_params[ch_index].norm(layers[c])
 
         # 2A) Image has 3 channels, no palette info, and no/only one cmap was given
         if palette is None and n_channels == 3 and not isinstance(render_params.cmap_params, list):
@@ -463,7 +646,7 @@ def _render_images(
                     "Consider using 'palette' instead."
                 )
 
-            _ax_show_and_transform(stacked, trans_data, ax, render_params.alpha)
+            _ax_show_and_transform(stacked, trans_data, ax, render_params.alpha, zorder=render_params.zorder)
 
         # 2B) Image has n channels, no palette/cmap info -> sample n categorical colors
         elif palette is None and not got_multiple_cmaps:
@@ -481,7 +664,7 @@ def _render_images(
             # Remove alpha channel so we can overwrite it from render_params.alpha
             colored = colored[:, :, :3]
 
-            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha)
+            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha, zorder=render_params.zorder)
 
         # 2C) Image has n channels and palette info
         elif palette is not None and not got_multiple_cmaps:
@@ -496,7 +679,7 @@ def _render_images(
             # Remove alpha channel so we can overwrite it from render_params.alpha
             colored = colored[:, :, :3]
 
-            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha)
+            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha, zorder=render_params.zorder)
 
         elif palette is None and got_multiple_cmaps:
             channel_cmaps = [cp.cmap for cp in render_params.cmap_params]  # type: ignore[union-attr]
@@ -509,7 +692,7 @@ def _render_images(
             # Remove alpha channel so we can overwrite it from render_params.alpha
             colored = colored[:, :, :3]
 
-            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha)
+            _ax_show_and_transform(colored, trans_data, ax, render_params.alpha, zorder=render_params.zorder)
 
         elif palette is not None and got_multiple_cmaps:
             raise ValueError("If 'palette' is provided, 'cmap' must be None.")
@@ -553,6 +736,7 @@ def _render_labels(
             scale=scale,
             is_label=True,
         )
+
     # rasterize spatial image if necessary to speed up performance
     if rasterize:
         label = _rasterize_if_necessary(
@@ -603,6 +787,17 @@ def _render_labels(
             seg_boundaries=render_params.outline,
             na_color=render_params.cmap_params.na_color,
         )
+        _cax = ax.imshow(
+            labels_infill,
+            rasterized=True,
+            cmap=None if categorical else render_params.cmap_params.cmap,
+            norm=None if categorical else render_params.cmap_params.norm,
+            alpha=render_params.fill_alpha,
+            origin="lower",
+            zorder=render_params.zorder,
+        )
+        _cax.set_transform(trans_data)
+        cax = ax.add_image(_cax)
 
         # Then overlay the contour
         labels_contour = _map_color_seg(
@@ -615,7 +810,6 @@ def _render_labels(
             seg_boundaries=render_params.outline,
             na_color=render_params.cmap_params.na_color,
         )
-
         _cax = ax.imshow(
             labels_contour,
             rasterized=True,
@@ -623,17 +817,9 @@ def _render_labels(
             norm=None if categorical else render_params.cmap_params.norm,
             alpha=render_params.outline_alpha,
             origin="lower",
+            zorder=render_params.zorder,
         )
-        _cax = ax.imshow(
-            labels_infill,
-            rasterized=True,
-            cmap=None if categorical else render_params.cmap_params.cmap,
-            norm=None if categorical else render_params.cmap_params.norm,
-            alpha=render_params.fill_alpha,
-            origin="lower",
-        )
-        _cax.set_transform(trans_data)
-        cax = ax.add_image(_cax)
+
     else:
         # Default: no alpha, contour = infill
         label = _map_color_seg(
@@ -654,6 +840,7 @@ def _render_labels(
             norm=None if categorical else render_params.cmap_params.norm,
             alpha=render_params.fill_alpha,
             origin="lower",
+            zorder=render_params.zorder,
         )
     _cax.set_transform(trans_data)
     cax = ax.add_image(_cax)
diff --git a/src/spatialdata_plot/pl/render_params.py b/src/spatialdata_plot/pl/render_params.py
index 2ccee43b..53b3a93e 100644
--- a/src/spatialdata_plot/pl/render_params.py
+++ b/src/spatialdata_plot/pl/render_params.py
@@ -80,6 +80,8 @@ class ShapesRenderParams:
     fill_alpha: float = 0.3
     scale: float = 1.0
     transfunc: Callable[[float], float] | None = None
+    method: str | None = None
+    zorder: int = 0
     table_name: str | None = None
 
 
@@ -96,6 +98,8 @@ class PointsRenderParams:
     alpha: float = 1.0
     size: float = 1.0
     transfunc: Callable[[float], float] | None = None
+    method: str | None = None
+    zorder: int = 0
     table_name: str | None = None
 
 
@@ -110,6 +114,7 @@ class ImageRenderParams:
     alpha: float = 1.0
     percentiles_for_norm: tuple[float | None, float | None] = (None, None)
     scale: str | None = None
+    zorder: int = 0
 
 
 @dataclass
@@ -128,3 +133,4 @@ class LabelsRenderParams:
     transfunc: Callable[[float], float] | None = None
     scale: str | None = None
     table_name: str | None = None
+    zorder: int = 0
diff --git a/src/spatialdata_plot/pl/utils.py b/src/spatialdata_plot/pl/utils.py
index 6fc3fa1a..6cc8f624 100644
--- a/src/spatialdata_plot/pl/utils.py
+++ b/src/spatialdata_plot/pl/utils.py
@@ -45,6 +45,7 @@
 from pandas.api.types import CategoricalDtype
 from scanpy import settings
 from scanpy.plotting._tools.scatterplots import _add_categorical_legend
+from scanpy.plotting._utils import add_colors_for_categorical_sample_annotation
 from scanpy.plotting.palettes import default_20, default_28, default_102
 from shapely.geometry import LineString, Polygon
 from skimage.color import label2rgb
@@ -835,8 +836,6 @@ def _get_palette(
 def _maybe_set_colors(
     source: AnnData, target: AnnData, key: str, palette: str | ListedColormap | Cycler | Sequence[Any] | None = None
 ) -> None:
-    from scanpy.plotting._utils import add_colors_for_categorical_sample_annotation
-
     color_key = f"{key}_colors"
     try:
         if palette is not None:
@@ -1626,6 +1625,9 @@ def _type_check_params(param_dict: dict[str, Any], element_type: str) -> dict[st
     if param_dict.get("table_name") and not isinstance(param_dict["table_name"], str):
         raise TypeError("Parameter 'table_name' must be a string .")
 
+    if param_dict.get("method") not in ["matplotlib", "datashader", None]:
+        raise ValueError("If specified, parameter 'method' must be either 'matplotlib' or 'datashader'.")
+
     return param_dict
 
 
@@ -1635,7 +1637,7 @@ def _validate_label_render_params(
     cmap: list[Colormap | str] | Colormap | str | None,
     color: str | None,
     fill_alpha: float | int,
-    contour_px: int,
+    contour_px: int | None,
     outline: bool,
     groups: list[str] | str | None,
     palette: list[str] | str | None,
@@ -1758,6 +1760,7 @@ def _validate_shape_render_params(
     norm: Normalize | None,
     scale: float | int,
     table_name: str | None,
+    method: str | None,
 ) -> dict[str, dict[str, Any]]:
     param_dict: dict[str, Any] = {
         "sdata": sdata,
@@ -1775,6 +1778,7 @@ def _validate_shape_render_params(
         "norm": norm,
         "scale": scale,
         "table_name": table_name,
+        "method": method,
     }
     param_dict = _type_check_params(param_dict, "shapes")
 
@@ -1804,6 +1808,7 @@ def _validate_shape_render_params(
 
         element_params[el]["palette"] = param_dict["palette"] if param_dict["col_for_color"] is not None else None
         element_params[el]["groups"] = param_dict["groups"] if param_dict["col_for_color"] is not None else None
+        element_params[el]["method"] = param_dict["method"]
 
     return element_params
 
@@ -1940,16 +1945,19 @@ def _ax_show_and_transform(
     ax: Axes,
     alpha: float | None = None,
     cmap: ListedColormap | LinearSegmentedColormap | None = None,
+    zorder: int = 0,
 ) -> None:
     if not cmap and alpha is not None:
         im = ax.imshow(
             array,
             alpha=alpha,
+            zorder=zorder,
         )
         im.set_transform(trans_data)
     else:
         im = ax.imshow(
             array,
             cmap=cmap,
+            zorder=zorder,
         )
         im.set_transform(trans_data)
diff --git a/tests/_images/Images_can_stick_to_zorder.png b/tests/_images/Images_can_stick_to_zorder.png
new file mode 100644
index 00000000..521e2bb9
Binary files /dev/null and b/tests/_images/Images_can_stick_to_zorder.png differ
diff --git a/tests/_images/Points_datashader_can_color_by_category.png b/tests/_images/Points_datashader_can_color_by_category.png
new file mode 100644
index 00000000..4a706d1f
Binary files /dev/null and b/tests/_images/Points_datashader_can_color_by_category.png differ
diff --git a/tests/_images/Points_datashader_continuous_color.png b/tests/_images/Points_datashader_continuous_color.png
new file mode 100644
index 00000000..cfd4e797
Binary files /dev/null and b/tests/_images/Points_datashader_continuous_color.png differ
diff --git a/tests/_images/Points_datashader_matplotlib_stack.png b/tests/_images/Points_datashader_matplotlib_stack.png
new file mode 100644
index 00000000..58fe248e
Binary files /dev/null and b/tests/_images/Points_datashader_matplotlib_stack.png differ
diff --git a/tests/_images/Shapes_datashader_can_color_by_category.png b/tests/_images/Shapes_datashader_can_color_by_category.png
new file mode 100644
index 00000000..fcba4711
Binary files /dev/null and b/tests/_images/Shapes_datashader_can_color_by_category.png differ
diff --git a/tests/_images/Shapes_datashader_can_color_by_value.png b/tests/_images/Shapes_datashader_can_color_by_value.png
new file mode 100644
index 00000000..5645d159
Binary files /dev/null and b/tests/_images/Shapes_datashader_can_color_by_value.png differ
diff --git a/tests/_images/Shapes_datashader_can_render_shapes.png b/tests/_images/Shapes_datashader_can_render_shapes.png
new file mode 100644
index 00000000..399a2d8e
Binary files /dev/null and b/tests/_images/Shapes_datashader_can_render_shapes.png differ
diff --git a/tests/pl/test_render_images.py b/tests/pl/test_render_images.py
index 9fa29893..313eda15 100644
--- a/tests/pl/test_render_images.py
+++ b/tests/pl/test_render_images.py
@@ -119,3 +119,6 @@ def test_plot_can_stack_render_images(self, sdata_blobs: SpatialData):
             .pl.render_images(element="blobs_image", channel=1, palette="blue", alpha=0.5)
             .pl.show()
         )
+
+    def test_plot_can_stick_to_zorder(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_shapes().pl.render_images().pl.show()
diff --git a/tests/pl/test_render_points.py b/tests/pl/test_render_points.py
index a7311a3e..b829392c 100644
--- a/tests/pl/test_render_points.py
+++ b/tests/pl/test_render_points.py
@@ -89,3 +89,18 @@ def test_plot_points_categorical_color(self, sdata_blobs: SpatialData):
 
         sdata_blobs["other_table"].obs["category"] = sdata_blobs["other_table"].obs["category"].astype("category")
         sdata_blobs.pl.render_points("blobs_points", color="category").pl.show()
+
+    def test_plot_datashader_continuous_color(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_points(
+            element="blobs_points", size=40, color="instance_id", alpha=0.6, method="datashader"
+        ).pl.show()
+
+    def test_plot_datashader_matplotlib_stack(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_points(
+            element="blobs_points", size=40, color="red", method="datashader"
+        ).pl.render_points(element="blobs_points", size=10, color="blue").pl.show()
+
+    def test_plot_datashader_can_color_by_category(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_points(
+            color="genes", groups="gene_b", palette="lightgreen", size=20, method="datashader"
+        ).pl.show()
diff --git a/tests/pl/test_render_shapes.py b/tests/pl/test_render_shapes.py
index 0d72bde4..328121df 100644
--- a/tests/pl/test_render_shapes.py
+++ b/tests/pl/test_render_shapes.py
@@ -286,3 +286,25 @@ def test_plot_shapes_categorical_color(self, sdata_blobs: SpatialData):
 
         sdata_blobs["table"].obs["category"] = sdata_blobs["table"].obs["category"].astype("category")
         sdata_blobs.pl.render_shapes("blobs_polygons", color="category").pl.show()
+
+    def test_plot_datashader_can_render_shapes(self, sdata_blobs: SpatialData):
+        sdata_blobs.pl.render_shapes(method="datashader").pl.show()
+
+    def test_plot_datashader_can_color_by_category(self, sdata_blobs: SpatialData):
+        RNG = np.random.default_rng(seed=42)
+        n_obs = len(sdata_blobs["blobs_polygons"])
+        adata = AnnData(RNG.normal(size=(n_obs, 10)))
+        adata.obs = pd.DataFrame(RNG.normal(size=(n_obs, 3)), columns=["a", "b", "c"])
+        adata.obs["category"] = RNG.choice(["a", "b", "c"], size=adata.n_obs)
+        adata.obs["instance_id"] = list(range(adata.n_obs))
+        adata.obs["region"] = "blobs_polygons"
+        table = TableModel.parse(adata=adata, region_key="region", instance_key="instance_id", region="blobs_polygons")
+        sdata_blobs["table"] = table
+
+        sdata_blobs.pl.render_shapes(element="blobs_polygons", color="category", method="datashader").pl.show()
+
+    def test_plot_datashader_can_color_by_value(self, sdata_blobs: SpatialData):
+        sdata_blobs["table"].obs["region"] = ["blobs_polygons"] * sdata_blobs["table"].n_obs
+        sdata_blobs["table"].uns["spatialdata_attrs"]["region"] = "blobs_polygons"
+        sdata_blobs.shapes["blobs_polygons"]["value"] = [1, 10, 1, 20, 1]
+        sdata_blobs.pl.render_shapes(element="blobs_polygons", color="value", method="datashader").pl.show()