基于 OpenLayers 的等值线生成与渲染

本文介绍如何在 Web 地图上从栅格瓦片数据中实时提取等值线并渲染。整个流程包含五个核心步骤：瓦片数据加载与解码 → 构建统一采样网格 → Marching Squares 提取 → 线段连接 → Chaikin 平滑。

---

渲染效果

1. 整体架构

┌─────────────┐     ┌──────────────┐     ┌──────────────────┐
│  JPEG 瓦片   │────▶│  像素解码 +   │────▶│  统一采样网格     │
│  (带 EXIF)   │     │  物理量还原   │     │  Float32Array    │
└─────────────┘     └──────────────┘     └────────┬─────────┘
                                                   │
                                                   ▼
┌─────────────┐     ┌──────────────┐     ┌──────────────────┐
│  OpenLayers  │◀────│  Chaikin 平滑 │◀────│  Marching Squares│
│  VectorLayer │     │  + 坐标转换   │     │  + 线段连接      │
└─────────────┘     └──────────────┘     └──────────────────┘

2. 瓦片数据加载与解码

我们的数据源是 JPEG 格式的瓦片图片，遵循标准的 {z}/{x}/{y} 瓦片编址。每张图片的 EXIF 元数据中存储了物理量的值域范围（uMin, uMax, vMin, vMax），像素的 R/G 通道分别编码了 U/V 两个分量的归一化值。

2.1 解码流程

JPEG 文件
  │
  ├── 解析 EXIF → 获取 { uMin, uMax, vMin, vMax }
  │
  └── 解码图像 → 256×256 RGBA 像素
        │
        └── 逐像素还原：
              u = uMin + (R / 255) × (uMax - uMin)
              v = vMin + (G / 255) × (vMax - vMin)
              speed = √(u² + v²)

2.2 陆地检测

对于海洋数据，需要区分陆地和海洋像素。判断条件：

const isLand = alpha < 128 || (red + green + blue) < 90

陆地像素的 speed 标记为 -1，后续 Marching Squares 会跳过包含陆地的网格单元。

2.3 高层级数据回退

数据通常只提供到某个最大缩放级别（如 maxZoom=2）。当地图缩放到更高级别时，需要回退到父级瓦片：

if (z > maxDataZoom) {
  const diff = z - maxDataZoom
  tz = maxDataZoom
  tx = Math.floor(x / 2^diff)
  ty = Math.floor(y / 2^diff)
}

3. 构建统一采样网格

等值线需要在连续的标量场上提取。如果对每个瓦片单独提取，瓦片边界处的等值线会断裂。因此我们将视口内所有瓦片的数据拼接成一个统一的大网格。

3.1 确定瓦片范围

根据当前视口的地图坐标范围，计算覆盖的瓦片索引：

const worldSize = 20037508.342789244 * 2  // EPSG:3857 全球宽度
const tileSize = worldSize / (1 << z)      // 每个瓦片的地图尺寸

const minTX = Math.floor((extent[0] - originX) / tileSize)
const maxTX = Math.floor((extent[2] - originX) / tileSize)
// ... 类似计算 Y 方向

注意 X 方向不做 clamp，这样可以支持 wrapX（跨日期变更线）。

3.2 动态采样步长

直接使用 256×256 的原始分辨率会导致网格过大。根据缩放级别动态调整采样步长：

缩放级别	步长 (step)	每瓦片采样点
≤ 2	4	64×64
≤ 4	2	128×128
> 4	1	256×256

同时设置网格维度上限（512×512），如果瓦片数量多导致超限，自动增大步长。

3.3 网格填充

const grid = new Float32Array(totalRows * totalCols).fill(-1)

// 对每个瓦片，将 speedData 按步长采样填入网格
for (let r = 0; r < sampledPerTile; r++) {
  for (let c = 0; c < sampledPerTile; c++) {
    const srcIdx = r * step * 256 + c * step
    grid[(gridRow + r) * totalCols + (gridCol + c)] = data.speedData[srcIdx]
  }
}

4. Marching Squares 算法

Marching Squares 是经典的二维等值线提取算法。它遍历网格中每个 2×2 的单元格，根据四个角点与阈值的关系确定等值线穿过该单元格的方式。

4.1 基本原理

对于一个 2×2 单元格，四个角点各有两种状态（≥ 阈值 或 < 阈值），共 2⁴ = 16 种组合：

角点编号与位权：

  TL(8) ─── TR(4)
   │          │
  BL(1) ─── BR(2)

caseIndex = (TL≥t ? 8:0) | (TR≥t ? 4:0) | (BR≥t ? 2:0) | (BL≥t ? 1:0)

4.2 16 种情况

case 0, 15: 无等值线（全在同侧）

case 1:  ╲     case 2:  ╱     case 3:  ──
         BL→L        B→R              L→R

case 4:  ╲     case 5:  ╳     case 6:  │
         T→R        鞍点             T→B

case 7:  ╱     case 8:  ╲     case 9:  │
         L→T        T→L              T→B

case 10: ╳     case 11: ╲     case 12: ──
         鞍点        T→R              L→R

case 13: ╱     case 14: ╲
         B→R        L→B

其中 case 5 和 case 10 是鞍点情况，一个单元格内有两条线段。

4.3 线性插值

等值线与单元格边的交点通过线性插值确定：

// 上边：TL 到 TR 之间
t = (threshold - TL) / (TR - TL)
交点 = (col + t, row)

// 右边：TR 到 BR 之间
t = (threshold - TR) / (BR - TR)
交点 = (col + 1, row + t)

4.4 性能优化

原始实现中每个线段创建一个 [x1, y1, x2, y2] 数组，大量小对象会给 GC 带来压力。优化方案：

// 使用 Float64Array 预分配连续内存
let buf = new Float64Array(capacity)
let count = 0

function push4(a, b, c, d) {
  const idx = count * 4
  if (idx + 4 > buf.length) {
    // 容量不足时翻倍扩容
    const nb = new Float64Array(buf.length * 2)
    nb.set(buf)
    buf = nb
  }
  buf[idx] = a; buf[idx+1] = b; buf[idx+2] = c; buf[idx+3] = d
  count++
}

同时内联 val() 和 lerp() 函数调用，预计算行偏移量，减少重复计算。

5. 线段连接

Marching Squares 输出的是大量独立的短线段。要得到连续的等值线，需要将共享端点的线段连接起来。

5.1 空间哈希索引

为了高效查找共享端点的线段，建立端点哈希索引：

function hashPt(x, y) {
  return `${Math.round(x * 1000)},${Math.round(y * 1000)}`
}

// 为每个线段的两个端点建立索引
for (let i = 0; i < segCount; i++) {
  const h1 = hashPt(seg[i].x1, seg[i].y1)
  const h2 = hashPt(seg[i].x2, seg[i].y2)
  index.get(h1).push(i)
  index.get(h2).push(i)
}

坐标乘以 1000 后取整作为哈希键，既能容忍浮点误差，又能快速匹配。

5.2 双向延伸

从一个未使用的线段出发，分别向头部和尾部延伸：

初始线段: A ──── B

向尾部延伸: 找到端点 B 匹配的下一个线段
  A ──── B ──── C ──── D

向头部延伸: 找到端点 A 匹配的线段
  F ──── E ──── A ──── B ──── C ──── D

5.3 避免 O(n²) 的 unshift

向头部延伸时，如果用 array.unshift() 插入，每次都是 O(n) 操作。优化方案是用一个单独的 headParts 数组收集头部点，最后一次性 reverse() 再拼接：

const headParts = []  // 收集头部延伸的点
const tailParts = []  // 收集尾部延伸的点

// 头部延伸时 push 到 headParts（反向）
headParts.push(newPoint)

// 最后拼接
headParts.reverse()
const line = headParts.concat(tailParts)

6. Chaikin 平滑

Marching Squares 产生的等值线是折线，视觉上不够平滑。Chaikin 细分算法通过迭代地在每条线段的 1/4 和 3/4 处插入新点来实现平滑。

6.1 算法原理

原始点:  P0 ──────── P1 ──────── P2

一次迭代后:
  P0 ── Q0 ── R0 ── Q1 ── R1 ── P2

其中:
  Q0 = 0.75 × P0 + 0.25 × P1  (P0→P1 的 1/4 处)
  R0 = 0.25 × P0 + 0.75 × P1  (P0→P1 的 3/4 处)
  Q1 = 0.75 × P1 + 0.25 × P2
  R1 = 0.25 × P1 + 0.75 × P2

6.2 实现

function chaikinSmooth(coords, iterations) {
  if (iterations <= 0 || coords.length < 3) return coords
  let pts = coords
  for (let iter = 0; iter < iterations; iter++) {
    const next = [pts[0]]  // 保留起点
    for (let i = 0; i < pts.length - 1; i++) {
      const p0 = pts[i], p1 = pts[i + 1]
      next.push([p0[0]*0.75 + p1[0]*0.25, p0[1]*0.75 + p1[1]*0.25])
      next.push([p0[0]*0.25 + p1[0]*0.75, p0[1]*0.25 + p1[1]*0.75])
    }
    next.push(pts[pts.length - 1])  // 保留终点
    pts = next
  }
  return pts
}

6.3 迭代次数的影响

迭代次数	点数增长	效果
0	1×	原始折线
1	~2×	略微圆滑
3	~8×	自然平滑（推荐）
5	~32×	非常平滑，但点数多
8	~256×	过度平滑，性能下降

7. 坐标转换与渲染

网格坐标（行列号）需要转换为地图坐标（EPSG:3857 米）：

// 网格左上角的地图坐标
const gridLeft = originX + minTX * tileSize
const gridTop  = originY - minTY * tileSize

// 每个网格单元的地图尺寸
const cellSize = (tileSize / 256) * step

// 网格坐标 → 地图坐标
mapX = gridLeft + col * cellSize
mapY = gridTop  - row * cellSize  // Y 轴方向相反

转换后的坐标用于创建 OpenLayers 的 LineString 几何体，添加到 VectorSource 中渲染。

8. 性能优化策略

8.1 结果缓存

以瓦片范围 + 步长 + 平滑度作为缓存键，视口不变时直接返回缓存结果：

const cacheKey = `${minTX},${maxTX},${minTY},${maxTY},${step},${smoothLevel}`
if (cacheKey === lastKey) return cachedFeatures

8.2 版本号取消

快速连续缩放时，旧的计算任务应该被取消：

async function updateContours() {
  const ver = ++version
  const features = await generate(ver)
  // 如果在计算期间又触发了新的更新，丢弃当前结果
  if (ver !== version) return
  source.clear()
  source.addFeatures(features)
}

8.3 样式对象复用

10 个等值线级别只创建 10 个 Style 对象，所有同级别的 Feature 共享同一个 Style 引用，避免重复创建。

8.4 Debounce

地图 moveend 事件触发后延迟 500ms 再执行更新，避免缩放过程中频繁计算。

9. 完整调用示例

import { useContourLayer } from '@/components/contourLayer'

// 在地图初始化后调用
const contour = useContourLayer(map, {
  url: '/data/{z}/{x}/{y}/current.jpg',
  maxDataZoom: 2,
  levels: [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0],
  colors: {
    0.1: '#0000ff', 0.2: '#0066ff', 0.3: '#00ccff',
    0.4: '#00ffcc', 0.5: '#00ff66', 0.6: '#66ff00',
    0.7: '#ccff00', 0.8: '#ffcc00', 0.9: '#ff6600',
    1.0: '#ff0000',
  },
  debounce: 500,
})

// 添加矢量图层
map.addLayer(new VectorLayer({
  source: contour.source,
  zIndex: 10,
}))

// 调整平滑度
contour.smoothLevel.value = 3

// 色斑图可以共用 tileLoader
const loader = contour.tileLoader

// 销毁
contour.destroy()

10. 总结

步骤	核心算法	输入	输出
数据加载	EXIF 解析 + 像素还原	JPEG 瓦片	Float32Array (speed)
网格构建	多瓦片拼接 + 动态采样	多个瓦片数据	统一 Float32Array 网格
等值线提取	Marching Squares	网格 + 阈值	线段集合 (Float64Array)
线段连接	空间哈希 + 双向延伸	线段集合	连续折线数组
平滑	Chaikin 细分	折线坐标	平滑曲线坐标
渲染	OpenLayers VectorLayer	Feature[]	地图上的等值线

整个流程的关键在于：将多个瓦片拼成统一网格保证跨瓦片连续性，Marching Squares + 线段连接从标量场中提取拓扑正确的等值线，Chaikin 平滑提升视觉效果，以及缓存 + 版本取消 + debounce保证交互流畅。

需要源码请联系