density contours stroke & fill (#948)

* density contours as an initializer: first pass

* density by z, carrying styles and channels (no reducer, only "first")

* fill: "density"

* density weight

* consistent thresholds across facets and series

note: we have to go around this bug in d3-contour: d3/d3-contour#57

* allow initializer composition; move initializer to the class; clean up

* density weight example

* error if x or y is undefined

(in the future we might route those to 1-dimensional transforms—KDE)

* document

* reduce img

* * don't apply the scales if we are already in pixel space (e.g. when composing with the hexbin transform)
* avoids a crash when there is no contour

* 1d density contours with frameAnchor

* adopt d3@7.4.5 for https://github.com/d3/d3-contour/releases/tag/v3.0.2

(not entirely sure if it's good practice to have all the yarn.lock changes in the PR, beyond those that are relevant)

* replace example image

* mike’s edits

* fix image dimensions

* isDensity

* distinct

* tweak tests

Co-authored-by: Mike Bostock <mbostock@gmail.com>

Author: Fil

README.md

@@ -1022,6 +1022,22 @@ Draws a mesh for the cell boundaries of the Voronoi tesselation of the points gi
 
 If a **z** channel is specified, the input points are grouped by *z*, and separate Voronoi tesselations are constructed for each group.
 
+### Density
+
+[<img src="./img/density-contours.png" width="320" height="200" alt="A scatterplot showing the relationship between the idle duration and eruption duration for Old Faithful">](https://observablehq.com/@observablehq/plot-density)
+
+[Source](./src/marks/density.js) · [Examples](https://observablehq.com/@observablehq/plot-density) · Draws regions of a two-dimensional point distribution in which the number of points per unit of screen space exceeds a certain density.
+
+#### Plot.density(*data*, *options*)
+
+Draws a region for each density level where the number of points given by the **x** and **y** channels, and possibly weighted by the **weight** channel, exceeds the given level. The **thresholds** option, which defaults to 20, indicates the approximate number of levels that will be computed at even intervals between 0 and the maximum density.
+
+If a **z**, **stroke** or **fill** channel is specified, the input points are grouped by series, and separate sets of contours are generated for each series.
+
+If stroke or fill is specified as *density*, a color channel is returned with values representing the density normalized between 0 and 1.
+
+If either of the **x** or **y** channels are not specified, the corresponding position is controlled by the **frameAnchor** option.
+
 ### Dot
 
 [<img src="./img/dot.png" width="320" height="198" alt="a scatterplot">](https://observablehq.com/@observablehq/plot-dot)
@@ -1135,7 +1151,7 @@ Returns a new image with the given *data* and *options*. If neither the **x** no
 
 ### Linear regression
 
-[<img src="./img/linear-regression.png" width="600" alt="a scatterplot of penguin culmens, showing the length and depth of several species, with linear regression models by species and for the whole population, illustrating Simpson’s paradox">](https://observablehq.com/@observablehq/plot-linear-regression)
+[<img src="./img/linear-regression.png" width="320" height="200" alt="a scatterplot of penguin culmens, showing the length and depth of several species, with linear regression models by species and for the whole population, illustrating Simpson’s paradox">](https://observablehq.com/@observablehq/plot-linear-regression)
 
 [Source](./src/marks/linearRegression.js) · [Examples](https://observablehq.com/@observablehq/plot-linear-regression) · Draws [linear regression](https://en.wikipedia.org/wiki/Linear_regression) lines with confidence bands, representing the estimated relation of a dependent variable (typically *y*) on an independent variable (typically *x*). The linear regression line is fit using the [least squares](https://en.wikipedia.org/wiki/Least_squares) approach. See Torben Jansen’s [“Linear regression with confidence bands”](https://observablehq.com/@toja/linear-regression-with-confidence-bands) and [this StatExchange question](https://stats.stackexchange.com/questions/101318/understanding-shape-and-calculation-of-confidence-bands-in-linear-regression) for details on the confidence interval calculation.
 

img/density-contours.png

@@ -50,7 +50,7 @@
     "vite": "2"
   },
   "dependencies": {
-    "d3": "^7.3.0",
+    "d3": "^7.4.5",
     "interval-tree-1d": "1",
     "isoformat": "0.2"
   },

img/linear-regression.png

@@ -1,7 +1,9 @@
 import {contourDensity, create, geoPath} from "d3";
-import {constant, maybeTuple} from "../options.js";
+import {identity, maybeTuple, maybeZ, valueof} from "../options.js";
 import {Mark} from "../plot.js";
-import {applyFrameAnchor, applyIndirectStyles, applyTransform} from "../style.js";
+import {coerceNumbers} from "../scales.js";
+import {applyFrameAnchor, applyDirectStyles, applyIndirectStyles, applyChannelStyles, applyTransform, distinct, groupZ} from "../style.js";
+import {initializer} from "../transforms/basic.js";
 
 const defaults = {
   ariaLabel: "density",
@@ -11,44 +13,151 @@ const defaults = {
 };
 
 export class Density extends Mark {
-  constructor(data, options = {}) {
-    const {x, y, bandwidth = 20, thresholds = 20} = options;
+  constructor(data, {x, y, z, weight, fill, stroke, ...options} = {}) {
+    // If fill or stroke is specified as “density”, then temporarily treat these
+    // as a literal color when computing defaults and maybeZ; below, we’ll unset
+    // these constant colors back to undefined since they will instead be
+    // populated by a channel generated by the initializer.
+    const fillDensity = isDensity(fill) && (fill = "currentColor", true);
+    const strokeDensity = isDensity(stroke) && (stroke = "currentColor", true);
     super(
       data,
       [
         {name: "x", value: x, scale: "x", optional: true},
-        {name: "y", value: y, scale: "y", optional: true}
+        {name: "y", value: y, scale: "y", optional: true},
+        {name: "z", value: maybeZ({z, fill, stroke}), optional: true},
+        {name: "weight", value: weight, optional: true}
       ],
-      options,
+      densityInitializer({...options, fill, stroke}, fillDensity, strokeDensity),
       defaults
     );
-    this.bandwidth = +bandwidth;
-    this.thresholds = +thresholds;
+    if (fillDensity) this.fill = undefined;
+    if (strokeDensity) this.stroke = undefined;
+    this.z = z;
+  }
+  filter(index) {
+    return index; // don’t filter contours constructed by initializer
   }
   render(index, scales, channels, dimensions) {
-    const {x: X, y: Y} = channels;
-    const {bandwidth, thresholds} = this;
-    const [cx, cy] = applyFrameAnchor(this, dimensions);
-    const {width, height} = dimensions;
+    const {contours} = channels;
+    const path = geoPath();
     return create("svg:g")
         .call(applyIndirectStyles, this, scales, dimensions)
         .call(applyTransform, this, scales)
-        .call(g => g.selectAll("path")
-          .data(contourDensity()
-              .x(X ? i => X[i] : constant(cx))
-              .y(Y ? i => Y[i] : constant(cy))
-              .size([width, height])
-              .bandwidth(bandwidth)
-              .thresholds(thresholds)
-            (index))
+        .call(g => g.selectAll()
+          .data(index)
           .enter()
           .append("path")
-            .attr("d", geoPath()))
-      .node();
+            .call(applyDirectStyles, this)
+            .call(applyChannelStyles, this, channels)
+            .attr("d", i => path(contours[i])))
+        .node();
   }
 }
 
 export function density(data, {x, y, ...options} = {}) {
   ([x, y] = maybeTuple(x, y));
   return new Density(data, {...options, x, y});
 }
+
+const dropChannels = new Set(["x", "y", "z", "weight"]);
+
+function densityInitializer(options, fillDensity, strokeDensity) {
+  let {bandwidth, thresholds} = options;
+  bandwidth = bandwidth === undefined ? 20 : +bandwidth;
+  thresholds = thresholds === undefined ? 20 : +thresholds; // TODO Allow an array of thresholds?
+  return initializer(options, function(data, facets, channels, scales, dimensions) {
+    const X = channels.x ? coerceNumbers(valueof(channels.x.value, scales[channels.x.scale] || identity)) : null;
+    const Y = channels.y ? coerceNumbers(valueof(channels.y.value, scales[channels.y.scale] || identity)) : null;
+    const W = channels.weight ? coerceNumbers(channels.weight.value) : null;
+    const Z = channels.z?.value;
+    const {z} = this;
+    const [cx, cy] = applyFrameAnchor(this, dimensions);
+    const {width, height} = dimensions;
+
+    // Group any of the input channels according to the first index associated
+    // with each z-series or facet. Drop any channels not be needed for
+    // rendering after the contours are computed.
+    const newChannels = Object.fromEntries(Object.entries(channels)
+      .filter(([key]) => !dropChannels.has(key))
+      .map(([key, channel]) => [key, {...channel, value: []}]));
+
+    // If the fill or stroke encodes density, construct new output channels.
+    const FD = fillDensity && [];
+    const SD = strokeDensity && [];
+    const k = 100; // arbitrary scale factor for readability
+
+    const density = contourDensity()
+        .x(X ? i => X[i] : cx)
+        .y(Y ? i => Y[i] : cy)
+        .weight(W ? i => W[i] : 1)
+        .size([width, height])
+        .bandwidth(bandwidth)
+        .thresholds(thresholds);
+
+    // If there are multiple facets or multiple series, first compute the
+    // contours for each facet-series independently; choose the set of contours
+    // with the maximum threshold value (density), and then apply this set’s
+    // thresholds to all the other facet-series. TODO With API changes to
+    // d3-contour, we could avoid recomputing the blurred grid and cache
+    // individual contours, making this more efficient.
+    if (facets.length > 1 || Z && facets.length > 0 && distinct(facets[0], Z)) {
+      let maxValue = 0;
+      let maxContours = [];
+      for (const facet of facets) {
+        for (const index of Z ? groupZ(facet, Z, z) : [facet]) {
+          const C = density(index);
+          if (C.length > 0) {
+            const c = C[C.length - 1];
+            if (c.value > maxValue) {
+              maxValue = c.value;
+              maxContours = C;
+            }
+          }
+        }
+      }
+      density.thresholds(maxContours.map(c => c.value));
+    }
+
+    // Generate contours for each facet-series.
+    const newFacets = [];
+    const contours = [];
+    for (const facet of facets) {
+      const contourFacet = [];
+      newFacets.push(contourFacet);
+      for (const index of Z ? groupZ(facet, Z, z) : [facet]) {
+        for (const contour of density(index)) {
+          contourFacet.push(contours.length);
+          contours.push(contour);
+          if (FD) FD.push(contour.value * k);
+          if (SD) SD.push(contour.value * k);
+          for (const key in newChannels) {
+            newChannels[key].value.push(channels[key].value[index[0]]);
+          }
+        }
+      }
+    }
+
+    // If the fill or stroke encodes density, ensure that a zero value is
+    // included so that the default color scale domain starts at zero. Otherwise
+    // if the starting range value is the same as the background color, the
+    // first contour might not be visible.
+    if (FD) FD.push(0);
+    if (SD) SD.push(0);
+
+    return {
+      data,
+      facets: newFacets,
+      channels: {
+        ...newChannels,
+        ...FD && {fill: {value: FD, scale: "color"}},
+        ...SD && {stroke: {value: SD, scale: "color"}},
+        contours: {value: contours}
+      }
+    };
+  });
+}
+
+function isDensity(value) {
+  return /^density$/i.test(value);
+}

package.json

@@ -184,6 +184,16 @@ function groupAesthetics({ariaLabel: AL, title: T, fill: F, fillOpacity: FO, str
   return [AL, T, F, FO, S, SO, SW, O, H].filter(c => c !== undefined);
 }
 
+export function distinct(I, X) {
+  const x = keyof(X[0]);
+  for (const i of I) {
+    if (keyof(X[i]) !== x) {
+      return true;
+    }
+  }
+  return false;
+}
+
 export function groupZ(I, Z, z) {
   const G = group(I, i => Z[i]);
   if (z === undefined && G.size > I.length >> 1) {