Past a few hundred thousand points, drawing one marker each stops working. Markers pile on top of markers, the file balloons, and a solid blob hides whatever structure the data had. Drawing the markers faster does not fix this, because the picture itself is wrong: overplotting throws away density.
datashade() takes the approach datashader made popular: do not draw
markers at all. Bin the points into a canvas-sized grid in one pass,
count how many land in each cell, and colour each cell by that count.
The cost is decoupled from both the number of points and the amount of
overplotting, and the result is honest about density.
datashade() returns a single raster_grob() you
draw like anything else.
The examples here are strange attractors, the same gallery shipped in
system.file("examples/attractors.R", package = "vellum").
They are a good stress test: a chaotic map iterated for millions of
points, with density spanning many orders of magnitude between the
sparse outer wisps and the bright folded core.
Generating an orbit
An attractor orbit is sequential: each point is a function of the one
before it, so it cannot be vectorised in R. vellum ships a small Rust
kernel for it, rs_attractor(), which returns 10 million
points in a fraction of a second. It is internal (hence the
:::), exposed here only to feed the example; the subject of
this article is what happens after the cloud exists.
n_points <- 2e6
attractor <- function(kind, p, n = n_points, x0 = 0.1, y0 = 0.1) {
n <- as.integer(n)
# jarl-ignore internal_function: the Rust kernel is internal; the vignette drives it directly
v <- vellum:::rs_attractor(kind, n, p[1], p[2], p[3], p[4], x0, y0)
list(x = v[seq_len(n)], y = v[n + seq_len(n)])
}
pts <- attractor("clifford", c(-1.4, 1.6, 1.0, 0.7))
length(pts$x)
#> [1] 2000000Shading in one call
Hand the cloud to datashade(). width and
height are the aggregation grid in cells, which is also the
output raster in pixels, and colors is the low-to-high
density ramp. A white first stop lets the sparse regions fade into the
page.
blues <- c("#ffffff", "#c6dbef", "#6baed6", "#2171b5", "#08306b")
vl_scene(5, 5, bg = "white") |>
draw(datashade(
pts$x, pts$y,
width = 500, height = 500,
colors = blues
))
Two million points, and the work scaled with the half-million grid cells, not the points. The same call handles twenty million with no change to the picture’s size on disk: still one 500x500 raster.
Why aggregation wins
The whole cost is one linear pass over the points to bin them, then a
colour lookup per grid cell. Concretely, datashade() calls
an aggregation kernel that is the O(N) heart of the method:
# jarl-ignore internal_function: the Rust kernel is internal; the vignette drives it directly
counts <- vellum:::rs_aggregate_2d(
pts$x, pts$y,
NULL, # optional per-point weights
120L, 120L, # grid size
min(pts$x), max(pts$x), min(pts$y), max(pts$y)
)
dim(counts)
#> NULL
range(counts) # the densest cell holds this many points
#> [1] 0 4985
mean(counts > 0) # fraction of the grid the orbit actually touches
#> [1] 0.5379167That count grid is all the downstream cost depends on. Whether it came from two million points or two hundred million, the shading step sees the same 120x120 matrix. Overplotting, the thing that ruins a scatter of markers, is exactly what the count measures here rather than hides.
Mapping density to colour
The densest cell of an attractor can hold thousands of times more
points than a faint outer cell. Map that range linearly and the bright
core saturates while everything faint collapses to nearly the
background. The how argument controls the density-to-colour
mapping; datashader’s default, "eq_hist" (histogram
equalisation), allocates colour by rank so structure stays visible
across the whole range.
ramp <- c("#ffffff", "#08306b")
panel <- function(how) {
datashade(pts$x, pts$y, width = 300, height = 300, colors = ramp, how = how)
}
vl_scene(6, 6, bg = "white") |>
push(vl_viewport(layout = grid_layout(
widths = vl_unit(c(1, 1), "null"),
heights = vl_unit(c(1, 1), "null")
))) |>
push(vl_viewport(row = 1, col = 1)) |> draw(panel("linear")) |> pop() |>
push(vl_viewport(row = 1, col = 2)) |> draw(panel("log")) |> pop() |>
push(vl_viewport(row = 2, col = 1)) |> draw(panel("cbrt")) |> pop() |>
push(vl_viewport(row = 2, col = 2)) |> draw(panel("eq_hist")) |> pop()
Clockwise from top-left: "linear" (only the core
survives), "log", "cbrt", and
"eq_hist" (the fullest picture). "log" and
"cbrt" are useful middle grounds when you want a mapping
with a fixed analytic meaning rather than one that depends on the data’s
rank distribution.
Lining up with data axes
datashade() bins over xlim by
ylim and returns a raster that fills its viewport
(npc 0..1). To place it against data axes,
draw it inside a vl_viewport() whose xscale /
yscale match the same limits, then draw axes, labels, or
reference lines in "native" units in that viewport. For
crisp bins, match width / height to the
viewport’s pixel size and keep the default
interpolate = FALSE.
For the attractors we centre a square window on the orbit so a square cell maps without distortion:
window <- function(x, y, pad = 1.05) {
xr <- range(x)
yr <- range(y)
half <- max(diff(xr), diff(yr)) / 2 * pad
list(
xlim = mean(xr) + c(-half, half),
ylim = mean(yr) + c(-half, half)
)
}
w <- window(pts$x, pts$y)
str(w)
#> List of 2
#> $ xlim: num [1:2] -1.5 2.06
#> $ ylim: num [1:2] -1.67 1.9Passing that window as xlim / ylim keeps
the aspect honest regardless of how the raw orbit range differs between
the two axes.
A gallery
Because each shaded attractor is just a grob, a gallery is a layout
of them. Push a grid_layout() of "null"
(flexible) tracks and drop one datashade() raster into each
cell. Square cells (from the square window above) keep every map
undistorted.
gallery <- list(
list(kind = "clifford", p = c(-1.4, 1.6, 1.0, 0.7),
pal = c("#ffffff", "#c6dbef", "#6baed6", "#2171b5", "#08306b")),
list(kind = "dejong", p = c( 1.4, -2.3, 2.4, -2.1),
pal = c("#ffffff", "#fdd0a2", "#fd8d3c", "#d94801", "#7f2704")),
list(kind = "svensson", p = c( 1.5, -1.8, 1.6, 0.9),
pal = c("#ffffff", "#c7e9c0", "#74c476", "#238b45", "#00441b")),
list(kind = "clifford", p = c(-1.8, -2.0, -0.5, -0.9),
pal = c("#ffffff", "#dadaeb", "#9e9ac8", "#6a51a3", "#3f007d"))
)
ncol <- 2
nrow <- 2
cell_px <- 300
s <- vl_scene(6, 6, dpi = 100, bg = "white") |>
push(vl_viewport(layout = grid_layout(
widths = vl_unit(rep(1, ncol), "null"),
heights = vl_unit(rep(1, nrow), "null")
)))
for (i in seq_along(gallery)) {
a <- gallery[[i]]
orbit <- attractor(a$kind, a$p)
w <- window(orbit$x, orbit$y)
img <- datashade(
orbit$x, orbit$y,
width = cell_px, height = cell_px,
xlim = w$xlim, ylim = w$ylim,
colors = a$pal, how = "eq_hist"
)
row <- (i - 1) %/% ncol + 1
col <- (i - 1) %% ncol + 1
s <- s |> push(vl_viewport(row = row, col = col)) |> draw(img) |> pop()
}
s
The shipped example (inst/examples/attractors.R) goes
further: it draws twelve panels with random parameters, keeping
only orbits that fill a fair fraction of the canvas (rejection sampling,
since most random parameters collapse to a point or diverge), each under
a random colormap. Run it for a fresh gallery every time:
Lines: dense timeseries
Points are not the only thing that overplots. A few hundred stacked
timeseries, each with a few thousand samples, is already a million line
vertices — draw them as vector paths and you get a solid band that hides
where the traces actually concentrate. The same aggregate-then-shade
idea applies: rasterise the lines into a grid where each cell
accumulates the coverage of the lines crossing it, then colour by that
density. datashade_lines() does this.
It takes a connected polyline. Pass group to pack many
series into one call — the line breaks wherever the group id changes (an
NA coordinate breaks it too), so a whole panel of traces is
a single call and a single raster.
set.seed(1)
k <- 600 # series
m <- 800 # samples each
steps <- matrix(rnorm(k * m, sd = 0.4), m, k)
walks <- apply(steps, 2, cumsum) # k random walks, m samples
t <- rep(seq_len(m), times = k) # x: shared time axis
y <- as.vector(walks) # y: all traces concatenated
grp <- rep(seq_len(k), each = m) # which series each vertex belongs to
vl_scene(6, 4, bg = "white") |>
push(vl_viewport(xscale = range(t), yscale = range(y))) |>
draw(datashade_lines(
t, y, group = grp,
width = 600, height = 400,
colors = c("#ffffff", "#fdd0a2", "#fd8d3c", "#d94801", "#7f2704"),
how = "eq_hist"
))
The bright core is where the walks pile up (near the start, before they diffuse apart); the faint outer envelope is the handful of traces that wandered far. That structure is exactly what an opaque band of 480,000 vertices would have hidden. Coverage is anti-aliased and summed, so overlapping traces brighten honestly rather than saturating to a single flat colour.
Segments: network edges
datashade_segments() is the same engine for
independent segments — each element is one
(x0, y0) -> (x1, y1) line, not a connected path. The
motivating case is the edges of a large graph: a hairball of tens of
thousands of edges drawn as vectors is an unreadable blob, but
aggregated it becomes an edge-density field showing which
regions of the layout the edges actually run through.
set.seed(2)
# A toy "layout": two loose clusters of nodes with many edges between them.
n_nodes <- 400
nx <- c(rnorm(n_nodes / 2, -1), rnorm(n_nodes / 2, 1))
ny <- rnorm(n_nodes)
n_edges <- 8000
a <- sample(n_nodes, n_edges, replace = TRUE)
b <- sample(n_nodes, n_edges, replace = TRUE)
vl_scene(6, 4, bg = "white") |>
push(vl_viewport(xscale = range(nx), yscale = range(ny))) |>
draw(datashade_segments(
nx[a], ny[a], nx[b], ny[b],
width = 600, height = 400,
colors = c("#ffffff", "#c6dbef", "#6baed6", "#2171b5", "#08306b"),
how = "eq_hist"
))
Per-segment weight works like the point path’s per-point
weight: pass edge weights and cells accumulate summed weight instead of
a plain crossing count.
Keeping thin marks visible: spread
Aggregated lines are often a single pixel wide, and a lone pixel
disappears when the image is shown small or over a busy background.
spread() dilates each non-empty pixel over a small
neighbourhood — the aggregation is untouched, the marks just get a
little bolder. dynspread() chooses the radius automatically
from how dense the image already is (dense regions spread less), which
is what the spread = "auto" shortcut on the
datashade_* functions calls.
sparse <- datashade_segments(
nx[a][1:1500], ny[a][1:1500], nx[b][1:1500], ny[b][1:1500],
width = 300, height = 300, colors = c("#ffffff", "#08306b")
)
vl_scene(6, 3, bg = "white") |>
push(vl_viewport(layout = grid_layout(
widths = vl_unit(c(1, 1), "null"),
heights = vl_unit(1, "null")
))) |>
push(vl_viewport(row = 1, col = 1)) |> draw(sparse) |> pop() |>
push(vl_viewport(row = 1, col = 2)) |> draw(dynspread(sparse)) |> pop()
Left, the raw one-pixel edges; right, the same raster after
dynspread(). Pass spread = 2 (a fixed radius)
or spread = "auto" (dynspread) straight to
datashade_lines() / datashade_segments() to
fold this into the shading call.
Recap
-
datashade()bins a point cloud into a grid and colours cells by density, so cost tracks the grid size, not the point count, and overplotting becomes signal rather than noise. -
datashade_lines()anddatashade_segments()do the same for dense timeseries and network edges: an anti-aliased line rasteriser accumulates coverage per cell, so overlapping lines add rather than saturate. -
how = "eq_hist"keeps structure visible across orders of magnitude;"log","cbrt", and"linear"are simpler alternatives. -
spread()/dynspread()(orspread =on thedatashade_*calls) dilate thin rasterised marks so single-pixel lines stay visible. - The result is one
raster_grob(). Draw it inside a viewport whose scales matchxlim/ylimto align it with axes, or tile many of them in agrid_layout()for a gallery.
In a grammar
datashade() is the low-level engine: you supply the
coordinates, the limits, and the ramp. A grammar layer on top can wire
all of that up from a plot spec. vellumplot exposes it
as mark_datashade(), which bins straight from a data frame
and fits into a normal plot with scales, guides, and facets. Its datashading
article pushes the same engine to its limit, shading the full US
Census (about 306 million points) two ways. ```