contour.js

import {blur2, contours, geoPath, max, min, nice, range, ticks, thresholdSturges, scaleUtc} from "d3";
import {createChannels} from "../channel.js";
import {create} from "../context.js";
import {labelof, identity, arrayify, map, isTemporal} from "../options.js";
import {applyPosition} from "../projection.js";
import {applyChannelStyles, applyDirectStyles, applyIndirectStyles, applyTransform, styles} from "../style.js";
import {initializer} from "../transforms/basic.js";
import {maybeThresholds} from "../transforms/bin.js";
import {AbstractRaster, maybeTuples, rasterBounds, sampler} from "./raster.js";

const defaults = {
  ariaLabel: "contour",
  fill: "none",
  stroke: "currentColor",
  strokeMiterlimit: 1,
  pixelSize: 2
};

export class Contour extends AbstractRaster {
  constructor(data, {smooth = true, value, ...options} = {}) {
    const channels = styles({}, options, defaults);

    // If value is not specified explicitly, look for a channel to promote. If
    // more than one channel is present, throw an error. (To disambiguate,
    // specify the value option explicitly.)
    if (value === undefined) {
      for (const key in channels) {
        if (channels[key].value != null) {
          if (value !== undefined) throw new Error("ambiguous contour value");
          value = options[key];
          options[key] = "value";
        }
      }
    }

    // For any channel specified as the literal (contour threshold) "value"
    // (maybe because of the promotion above), propagate the label from the
    // original value definition.
    if (value != null) {
      const v = {transform: (D) => D.map((d) => d.value), label: labelof(value)};
      for (const key in channels) {
        if (options[key] === "value") {
          options[key] = v;
        }
      }
    }

    // If the data is null, then we’ll construct the raster grid by evaluating a
    // function for each point in a dense grid. The value channel is populated
    // by the sampler initializer, and hence is not passed to super to avoid
    // computing it before there’s data.
    if (data == null) {
      if (value == null) throw new Error("missing contour value");
      options = sampler("value", {value, ...options});
      value = null;
    }

    // Otherwise if data was provided, it represents a discrete set of spatial
    // samples (often a grid, but not necessarily). If no interpolation method
    // was specified, default to nearest.
    else {
      let {interpolate} = options;
      if (value === undefined) value = identity;
      if (interpolate === undefined) options.interpolate = "nearest";
    }

    // Wrap the options in our initializer that computes the contour geometries;
    // this runs after any other initializers (and transforms).
    super(data, {value: {value, optional: true}}, contourGeometry(options), defaults);

    // With the exception of the x, y, x1, y1, x2, y2, and value channels, this
    // mark’s channels are not evaluated on the initial data but rather on the
    // contour multipolygons generated in the initializer.
    const contourChannels = {geometry: {value: identity}};
    for (const key in this.channels) {
      const channel = this.channels[key];
      const {scale} = channel;
      if (scale === "x" || scale === "y" || key === "value") continue;
      contourChannels[key] = channel;
      delete this.channels[key];
    }
    this.contourChannels = contourChannels;
    this.smooth = !!smooth;
  }
  filter(index, {x, y, value, ...channels}, values) {
    // Only filter channels constructed by the contourGeometry initializer; the
    // x, y, and value channels must be filtered by the initializer itself.
    return super.filter(index, channels, values);
  }
  render(index, scales, channels, dimensions, context) {
    const {geometry: G} = channels;
    const path = geoPath();
    return create("svg:g", context)
      .call(applyIndirectStyles, this, dimensions, context)
      .call(applyTransform, this, scales)
      .call((g) => {
        g.selectAll()
          .data(index)
          .enter()
          .append("path")
          .call(applyDirectStyles, this)
          .attr("d", (i) => path(G[i]))
          .call(applyChannelStyles, this, channels);
      })
      .node();
  }
}

function contourGeometry({thresholds, interval, ...options}) {
  thresholds = maybeThresholds(thresholds, interval, thresholdSturges);
  return initializer(options, function (data, facets, channels, scales, dimensions, context) {
    const [x1, y1, x2, y2] = rasterBounds(channels, scales, dimensions, context);
    const dx = x2 - x1;
    const dy = y2 - y1;
    const {pixelSize: k, width: w = Math.round(Math.abs(dx) / k), height: h = Math.round(Math.abs(dy) / k)} = this;
    const kx = w / dx;
    const ky = h / dy;
    const temporal = isTemporal(channels.value.value);
    const V = temporal && this.blur > 0 ? Float64Array.from(channels.value.value) : channels.value.value;
    const VV = []; // V per facet

    // Interpolate the raster grid, as needed.
    if (this.interpolate) {
      const {x: X, y: Y} = applyPosition(channels, scales, context);
      // Convert scaled (screen) coordinates to grid (canvas) coordinates.
      const IX = map(X, (x) => (x - x1) * kx, Float64Array);
      const IY = map(Y, (y) => (y - y1) * ky, Float64Array);
      // The contour mark normally skips filtering on x, y, and value, so here
      // we’re careful to use different names (0, 1, 2) when filtering.
      const ichannels = [channels.x, channels.y, channels.value];
      const ivalues = [IX, IY, V];
      for (const facet of facets) {
        const index = this.filter(facet, ichannels, ivalues);
        VV.push(this.interpolate(index, w, h, IX, IY, V));
      }
    }

    // Otherwise, chop up the existing dense raster grid into facets, if needed.
    // V must be a dense grid in projected coordinates; if there are multiple
    // facets, then V must be laid out vertically as facet 0, 1, 2… etc.
    else if (facets) {
      const n = w * h;
      const m = facets.length;
      for (let i = 0; i < m; ++i) VV.push(V.slice(i * n, i * n + n));
    } else {
      VV.push(V);
    }

    // Blur the raster grid, if desired.
    if (this.blur > 0) for (const V of VV) blur2({data: V, width: w, height: h}, this.blur);

    // Compute the contour thresholds.
    const T = maybeTicks(thresholds, V, ...finiteExtent(VV), temporal);
    if (T === null) throw new Error(`unsupported thresholds: ${thresholds}`);

    // Compute the (maybe faceted) contours.
    const {contour} = contours().size([w, h]).smooth(this.smooth);
    const contourData = [];
    const contourFacets = [];
    for (const V of VV) {
      contourFacets.push(range(contourData.length, contourData.push(...map(T, (t) => contour(V, t)))));
    }

    // Rescale the contour multipolygon from grid to screen coordinates.
    for (const {coordinates} of contourData) {
      for (const rings of coordinates) {
        for (const ring of rings) {
          for (const point of ring) {
            point[0] = point[0] / kx + x1;
            point[1] = point[1] / ky + y1;
          }
        }
      }
    }

    // Compute the deferred channels.
    return {
      data: contourData,
      facets: contourFacets,
      channels: createChannels(this.contourChannels, contourData)
    };
  });
}

// Apply the thresholds interval, function, or count, and return an array of
// ticks. d3-contour unlike d3-array doesn’t pass the min and max automatically,
// so we do that here to normalize, and also so we can share consistent
// thresholds across facets. When an interval is used, note that the lowest
// threshold should be below (or equal) to the lowest value, or else some data
// will be missing.
function maybeTicks(thresholds, V, min, max, temporal) {
  if (typeof thresholds?.range === "function") return thresholds.range(thresholds.floor(min), max);
  if (typeof thresholds === "function") thresholds = thresholds(V, min, max);
  if (typeof thresholds !== "number") return arrayify(thresholds);
  if (temporal) return scaleUtc().domain([min, max]).nice(thresholds).ticks(thresholds);
  const tz = ticks(...nice(min, max, thresholds), thresholds);
  while (tz[tz.length - 1] >= max) tz.pop();
  while (tz[1] < min) tz.shift();
  return tz;
}

export function contour() {
  return new Contour(...maybeTuples("value", ...arguments));
}

function finiteExtent(VV) {
  return [min(VV, (V) => min(V, finite)), max(VV, (V) => max(V, finite))];
}

function finite(x) {
  return isFinite(x) ? x : NaN;
}