d3-contour 实战：用 marching squares 在矩形数值网格上计算等高线与核密度估计

本篇技术指南基于 d3 官方文档 docs/d3-contour.md 及其两篇子文档（Contours、Density estimation），完整覆盖 d3-contour 模块的两大核心能力：对矩形数值网格应用 marching squares 算法生成 GeoJSON MultiPolygon 等高线，以及基于快速二维核密度估计（KDE）的点云密度等高线。读完后你将掌握 contours() 与 contourDensity() 全部 API 的参数、取值与默认值，理解输入网格的坐标映射规则（values[i + jn] 与平面坐标 ⟨i + 0.5, j + 0.5⟩），并能结合 geoPath/geoIdentity 将等高线渲染到 SVG。

模块定位：把数值网格变成 GeoJSON 多边形

d3-contour 模块的核心作用是：对一个矩形数值网格应用 marching squares（移动方格）算法，计算等高线多边形。官方文档中的示例图展示了新西兰奥克兰的 Maungawhau（又称 Mt. Eden）火山的等高线地形，数据来自仓库内的真实数据集：

• 文档中的示例由 PlotRender 组件 驱动，数据加载逻辑见 volcano.data.js，它读取 volcano.json；
• 经核实该 JSON 结构为 { width: 87, height: 61, values: [...] }，values 长度 5307（即 87 × 61），与文档示例中传入 Plot.contour 的 width/height 选项完全对应；
• 在 d3 主包中，该模块通过 src/index.js 中的 export * from "d3-contour" 统一导出，依赖版本见 package.json：d3 7.9.0 依赖 d3-contour: ^4.0.2。

文档中的官方示例用法（Observable Plot 风格）如下，等价于先用 d3.contours() 生成多边形再传给图形库渲染：

Plot.contour(volcano.values, {
  width: volcano.width,   // 87
  height: volcano.height, // 61
  fill: Plot.identity,
  stroke: "black",
  interval: 5            // 阈值间隔
})

两个子主题的完整 API 文档分别在 Contours 子文档 和 Density estimation 子文档，下文完整继承并展开这两篇文档的全部内容。

等高线生成器 contours()

contours() 与 contours(values)

contours() 以默认设置构造一个新的等高线生成器：

const contours = d3.contours()
    .size([width, height])
    .thresholds([0, 1, 2, 3, 4]);

调用 contours(values) 时，对给定的 values 数组计算等高线，返回一个 GeoJSON MultiPolygon geometry 对象数组。每个 geometry 对象表示输入值大于或等于对应阈值值的区域；该阈值通过 geometry.value 暴露。

输入网格的坐标映射规则（这是正确使用本模块的关键，直接来自子文档）：

• values 必须是长度为 n × m 的数组，其中 [n, m] 是生成器的 size；
• 每个 values[i + jn] 表示位置 ⟨i, j⟩ 处的值；
• 多边形采用平面坐标（planar coordinates），网格元素 i + jn 对应坐标 ⟨i + 0.5, j + 0.5⟩。

子文档给出的一个完整可运行示例：对 Goldstein–Price 函数（全局优化测试函数）在 -2 ≤ x ≤ 2、-2 ≤ y ≤ 1 范围内采样构造 256×256 网格：

var n = 256, m = 256, values = new Array(n * m);
for (var j = 0.5, k = 0; j < m; ++j) {
  for (var i = 0.5; i < n; ++i, ++k) {
    values[k] = goldsteinPrice(i / n * 4 - 2, 1 - j / m * 3);
  }
}

function goldsteinPrice(x, y) {
  return (1 + Math.pow(x + y + 1, 2) * (19 - 14 * x + 3 * x * x - 14 * y + 6 * x * x + 3 * y * y))
      * (30 + Math.pow(2 * x - 3 * y, 2) * (18 - 32 * x + 12 * x * x + 48 * y - 36 * x * y + 27 * y * y));
}

注意循环中 i 和 j 都从 0.5 开始（而非 0），这正是网格中心点坐标 ⟨i + 0.5, j + 0.5⟩ 约定的体现。

渲染：交给 geoPath 与 geoIdentity

由于等高线多边形是标准 GeoJSON，可以借助常规工具进行变换与显示。典型做法是把返回的 geometry 传给 geoPath（参见 d3-geo/path.md），并将关联的投影设为 null 或 geoIdentity（平面坐标无需投影变换，geoIdentity 见 projection.md 中 ## geoIdentity() 一节）。

这意味着等高线结果还能进一步被投影变换与拼接（如文档提及的 geoProject、geoStitch，属于 d3-geo-projection 生态），例如把全球地表温度 GeoTIFF 的等高线显示在 Natural Earth 投影中。此外文档还列出两类典型应用场景：

1. 加载 GeoTIFF 地表温度生成全球等高线；
2. 对含噪声的灰度 PNG（如云量）做模糊后生成平滑等高线——模糊步骤可用 d3 生态中的 d3.blur 配合，网格数据即 0–1 灰度值。

生成器配置方法

contours.contour(values, threshold) — 计算单条等高线，返回一个 MultiPolygon geometry 对象，表示输入值大于或等于给定 threshold 的区域；geometry.value 即该阈值。输入网格的 values[i + jn] 映射规则与 contours(values) 相同。适合只关心某一条等值线（如"温度 ≥ 20°C 的区域"）的场景，避免计算全部阈值。

contours.size(size) — 设置输入 values 网格的预期尺寸。size 为数组 [n, m]：n 是列数，m 是行数，二者必须为正整数；不传参时返回当前 size，默认 [1, 1]。火山示例中即为 .size([87, 61])。

contours.smooth(smooth) — 设置是否对生成的等高线多边形做线性插值平滑。不传参时返回当前平滑标志，默认 true。关闭后多边形边会更贴近网格的阶梯状折线；开启后在相邻单元格值之间按比例内插交点，轮廓更光滑。

contours.thresholds(thresholds) — 设置阈值生成器（函数或数组），不传参时返回当前生成器，默认实现 Sturges 公式（thresholdSturges，参见 d3-array/bin.md）。阈值语义规则：

• 阈值定义为数组 [x0, x1, …]；
• 第一条等高线对应值 ≥ x0 的区域，第二条对应值 ≥ x1 的区域，依此类推；
• 因此每个阈值恰好生成一个 MultiPolygon geometry 对象，阈值通过 geometry.value 暴露；
• 若传入的是一个数量 count 而非数组，则对输入值的 extent 均匀划分为约 count 个区间生成阈值（基于 ticks 生成"好看"的刻度值）。

密度估计 contourDensity()

等高线还能展示点云的估计密度，这对避免大数据集上的过度绘图（overplotting）很有用——散点图在数万点时几乎糊成一片，而密度等高线能用嵌套的若干层多边形清晰表达分布形态。contourDensity 方法实现了快速的二维核密度估计：先用高斯核把散点"泼"到规则网格上得到密度值矩阵，再对该矩阵走与上文相同的 marching squares 流程。

典型示例：Old Faithful（老忠实间歇泉）空闲时长与喷发时长的关系散点，以及 53,940 颗钻石重量与价格的密度等高线——后者正是"点云密度代替散点"的动机所在。

contourDensity() 与 density(data)

contourDensity() 以默认设置构造新的密度估计器。调用 density(data) 对给定的 data 数组估计密度等高线，返回 GeoJSON MultiPolygon geometry 数组：

• 每个 geometry 表示估计的每平方像素点数大于或等于对应阈值值的区域，geometry.value 为该阈值；
• 返回的 geometry 通常传给 geoPath 显示，投影使用 null 或 geoIdentity；
• 每个数据点的 x/y 坐标由 density.x 与 density.y 计算，density.weight 指定每个点的相对贡献（默认 1）；
• 生成的等高线仅在估计器定义的 size 范围内准确——超出 size 边界的点会被截断处理。

估计器配置方法

density.x(x) — 设置 x 坐标访问器；不传参返回当前访问器，默认：

function x(d) {
  return d[0];
}

density.y(y) — 设置 y 坐标访问器；不传参返回当前访问器，默认：

function y(d) {
  return d[1];
}

即默认输入格式是 [x, y] 二元组数组；若数据是对象，需显式传入 d => d.latitude 之类的访问器。

density.weight(weight) — 设置点权重访问器，默认：

function weight() {
  return 1;
}

density.size(size) — 设置估计器的尺寸边界，size 为 [width, height]，其中 width 是最大 x 值、height 是最大 y 值；不传参返回当前 size，默认 [960, 500]。同样地，密度等高线只在定义的 size 范围内准确，因此 size 应覆盖数据的实际取值范围（例如钻石数据的 [maxCarat, maxPrice]）。

density.cellSize(cellSize) — 设置底层分箱网格中单个单元格的大小，必须为正整数；不传参返回当前值，默认 4。两个重要特性：

• 单元格大小会向下取整到最近的 2 的幂；
• 更小的单元格产生更精细的多边形，但计算代价更高——这是精度与性能的权衡旋钮。

density.thresholds(thresholds) — 设置阈值生成器（函数或数组），不传参返回当前生成器，默认生成约 20 个好看的密度阈值。语义与 contours.thresholds 一致：每个阈值恰好对应一个 MultiPolygon，geometry.value 为阈值；第一个阈值 x0 通常应大于零（密度值非负，从 0 开始会生成覆盖全域的无意义外环）。若传入数量 count 而非数组，则生成约 count 个均匀间隔的"好看"阈值（基于 ticks）。

density.bandwidth(bandwidth) — 设置高斯核的带宽（标准差）；不传参返回当前带宽，默认 20.4939…。指定值会被舍入到该实现当前支持的最接近值，且必须非负。带宽控制密度场的"平滑程度"：带宽越小密度估计越贴合局部点簇，越大越平滑。

density.contours(data) — 返回一个 contour(value) 函数，可在给定数据上计算任意单条密度等高线而无需重算底层网格；返回的 contour 函数还暴露 contour.max 值，表示网格中的最大密度。适合做交互式调整：预计算一次网格，然后按用户选择的密度值反复取单条等高线。

小结：从网格到图形的完整链路

综合两篇子文档，d3-contour 的使用路径可以归纳为一条固定链路：

1. 准备数据：要么直接拥有 n × m 的平面数值网格（按 values[i + jn] 排布），要么从点云出发用 contourDensity 先估计出密度网格；
2. 配置生成器：size 声明网格规模（contours）或数据域（density），thresholds 决定输出多边形层数（默认分别为 Sturges 公式与约 20 个"好看"阈值），smooth 控制插值平滑，cellSize/bandwidth 控制密度估计的分辨率与核宽；
3. 生成 GeoJSON：contours(values) / density(data) 输出 MultiPolygon 数组，每个对象带 value 字段标识阈值，contour(values, t) / density.contours(data) 输出单条等高线；
4. 渲染：把 geometry 交给 geoPath（投影取 null 或 geoIdentity）转成 SVG path，即可填充、描边或动画化。

由于输出是标准 GeoJSON，等高线还能无缝接入 d3-geo 的投影、变换与拼接工具链，从平面图表扩展到地理可视化。本文所有 API 细节、默认值与坐标约定均以仓库中 d3-contour 主文档、Contours 子文档 与 Density estimation 子文档 为准，适用前提为 d3 7.9.0（依赖 d3-contour ^4.0.2）。