Datashader doesn't speed up rendering of points

[timtreis]

To get a feeling for when rendering with datashader is faster than matplotlib, we prepared a notebook contrasting the two: https://github.com/scverse/spatialdata-notebooks/blob/main/notebooks/examples/speed_up_illustration.ipynb

However, in doing so we found that rendering points is not really sped up for some reason.

Why is that?

cc @Sonja-Stockhaus

[LucaMarconato]

Very strange, was this maybe due to the old implementation in which the figure size of datashader was given by the value of the coordinates/transformed coordinates? I think this has now been fixed. Also a way to test for that is to make random points with bound values (like points with coordinates in [0, 1]).

[timtreis]

Performance investigation findings

We profiled the datashader vs matplotlib rendering path for points and identified 5 bottlenecks that together made datashader consistently slower than matplotlib at all point counts.

Bottleneck 1: Dask DataFrame passed to datashader (~137x overhead)

PointsModel.parse() always returns a dask DataFrame, which is then passed to cvs.points(). Datashader's dask path triggers scheduler overhead on every aggregation call, even though the underlying data is already computed in memory.

# Warmed-up benchmark at 100K points:
cvs.points(pandas_df, "x", "y", agg=ds.count())   # 0.0005s
cvs.points(dask_df,   "x", "y", agg=ds.count())   # 0.0700s  (140x slower)

Fix: Call .compute() on the transformed element before datashader aggregation:

# render.py, datashader branch
transformed_element = PointsModel.parse(
    trans.transform(sdata_filt.points[element][["x", "y"]]),
    annotation=...,
    transformations={coordinate_system: Identity()},
).compute()  # <-- materialize to pandas

Bottleneck 2: Double extent computation (~0.05s at 500K)

The extent was computed via get_extent() on the dask DataFrame (triggering a .compute()), and then .compute() was called again for the aggregation. Both operations materialize the same data.

Fix: Compute the DataFrame once, then derive extent directly from the pandas DataFrame:

def _datashader_canvas_from_dataframe(df, ax, fig_params):
    x_ext = [min(0, float(df["x"].min())), float(df["x"].max())]
    y_ext = [min(0, float(df["y"].min())), float(df["y"].max())]
    return _compute_datashader_canvas_params(x_ext, y_ext, ax, fig_params)

Bottleneck 3: Per-point _hex_no_alpha() calls (~0.24s at 500K)

Line 909 in render.py called _hex_no_alpha() on every point in a list comprehension. In the common case (no alpha channel in hex colors), no stripping is needed at all.

# Before: O(n) string processing
color_vector = np.asarray([_hex_no_alpha(x) for x in color_vector])

# After: O(unique_colors) check, skip if not needed
unique_map = {c: _hex_no_alpha(c) for c in set(color_vector)}
if any(v != k for k, v in unique_map.items()):
    color_vector = np.asarray([unique_map[c] for c in color_vector])

Bottleneck 4: _build_datashader_color_key iterated all points (~0.13s at 100K)

This function iterated over every point to find the first color for each category, using expensive per-element pandas Index lookups (categories[code]). Additionally, the _inject_ds_nan_sentinel() adds a sentinel category that never appears in codes, preventing the early-exit condition from triggering.

Fix: Pre-extract as numpy arrays, count only present categories, early-exit:

categories_arr = np.asarray(cat_series.categories, dtype=str)
codes = np.asarray(cat_series.codes)
n_present = len(np.unique(codes[codes >= 0]))
first_color = {}
for code, color in zip(codes, color_vector, strict=False):
    if code < 0:
        continue
    cat_name = categories_arr[code]
    if cat_name not in first_color:
        first_color[cat_name] = _hex_no_alpha(color) if isinstance(color, str) and color.startswith("#") else color
        if len(first_color) == n_present:
            break

Bottleneck 5: _want_decorations O(n) set creation (~0.025s at 500K)

set(cv.tolist()) converted a 500K numpy array to a Python list then to a set. A numpy comparison is much faster:

# Before
unique_vals = set(cv.tolist())
if len(unique_vals) != 1:
    return True

# After
first = cv.flat[0]
if not (cv == first).all():
    return True

Results (after JIT warmup)

No coloring:

n	matplotlib	datashader (before)	datashader (after)	speedup
10K	0.065s	0.997s (0.07x)	0.055s	1.43x
50K	0.100s	0.207s (0.51x)	0.093s	1.20x
100K	0.145s	0.267s (0.58x)	0.111s	1.34x
500K	0.540s	0.763s (0.75x)	0.410s	1.37x
1M	1.040s	0.861s (1.21x)	0.780s	1.40x

Categorical coloring:

n	matplotlib	datashader (after)	speedup
100K	0.168s	0.154s	1.09x
500K	0.605s	0.381s	1.59x

At small sizes (<50K with categorical coloring), datashader's fixed initialization overhead (canvas setup, shading, image conversion) still dominates. This is inherent to datashader. The auto-selection threshold of 10K points remains appropriate.

Remaining shared overhead

Both paths share a common preprocessing pipeline (~0.3-0.4s at 500K) that includes dask compute, _reparse_points (pandas→dask→PointsModel round-trips), AnnData construction, and color resolution. Optimizing this shared path would benefit both backends but is a separate concern.