站内搜索

输入关键词开始搜索

从 Cesium 地形采样生成局部等高线

文章目录

Cesium 可以提供指定位置的地形高程,但不会直接返回任意区域的等高线。局部等高线需要把离散高程重新组织成规则网格,再从网格中提取相同高度的线。

完整管线是:

用户绘制多边形
  → 计算外接矩形
  → 构造规则采样网格
  → 批量采样地形高度
  → 平滑网格
  → Marching Squares 提取线段
  → 连接和裁剪线段
  → 批量渲染与标注

从绘制区域得到采样范围

先去掉重复闭合点,再计算经纬度范围:

function normalizePolygon(positions: Cesium.Cartesian3[]) {
  if (positions.length < 3) return []

  const first = positions[0]
  const last = positions[positions.length - 1]

  return Cesium.Cartesian3.distance(first, last) < 0.01
    ? positions.slice(0, -1)
    : positions
}

function calculateBounds(positions: Cesium.Cartesian3[]) {
  let west = Infinity
  let south = Infinity
  let east = -Infinity
  let north = -Infinity

  positions.forEach(position => {
    const point = Cesium.Cartographic.fromCartesian(position)
    const longitude = Cesium.Math.toDegrees(point.longitude)
    const latitude = Cesium.Math.toDegrees(point.latitude)

    west = Math.min(west, longitude)
    east = Math.max(east, longitude)
    south = Math.min(south, latitude)
    north = Math.max(north, latitude)
  })

  return { west, south, east, north }
}

采样范围比用户多边形外接矩形扩大约 10%。这样边界附近的等高线不会因为缺少外围网格而提前中断:

const paddingX = (bounds.east - bounds.west) * 0.1
const paddingY = (bounds.north - bounds.south) * 0.1

const expandedBounds = {
  west: bounds.west - paddingX,
  east: bounds.east + paddingX,
  south: bounds.south - paddingY,
  north: bounds.north + paddingY,
}

最终生成的线仍会裁剪回用户多边形,扩大的范围只用于计算。

构造规则高程网格

假设密度为 density,构造 density × density 个采样点:

function createSamplingGrid(bounds: RegionBounds, density: number) {
  const positions: Cesium.Cartographic[] = []

  for (let row = 0; row < density; row += 1) {
    for (let column = 0; column < density; column += 1) {
      const longitude = Cesium.Math.lerp(
        bounds.west,
        bounds.east,
        column / (density - 1),
      )
      const latitude = Cesium.Math.lerp(
        bounds.south,
        bounds.north,
        row / (density - 1),
      )

      positions.push(Cesium.Cartographic.fromDegrees(longitude, latitude))
    }
  }

  return positions
}

density = 200 会生成 40000 个采样点。密度增加一倍,点数会增加到四倍,因此 UI 不应把它当成普通线性参数。

批量采样地形高度

const positions = createSamplingGrid(expandedBounds, density)

const sampled = await Cesium.sampleTerrainMostDetailed(
  viewer.terrainProvider,
  positions,
)

把一维结果恢复成二维数组:

const data: number[][] = []
let index = 0

for (let row = 0; row < density; row += 1) {
  const values: number[] = []

  for (let column = 0; column < density; column += 1) {
    const sample = sampled[index]
    values.push(Number.isFinite(sample.height) ? sample.height : 0)
    index += 1
  }

  data.push(values)
}

生产代码还要处理请求失败和长时间无响应。超时只表示当前分析终止,不能真正取消已经发出的地形请求:

const result = await Promise.race([
  Cesium.sampleTerrainMostDetailed(viewer.terrainProvider, positions),
  new Promise<never>((_, reject) => {
    setTimeout(() => reject(new Error('terrain sample timeout')), 30000)
  }),
])

先平滑,再提取

地形采样中的局部噪声会产生很多短小碎线。可以使用一个 3×3 高斯核做一次平滑:

function smoothGrid(data: number[][]) {
  const rows = data.length
  const columns = data[0].length
  const result = data.map(row => row.slice())

  for (let i = 1; i < rows - 1; i += 1) {
    for (let j = 1; j < columns - 1; j += 1) {
      const center = data[i][j]
      const north = data[i - 1][j]
      const south = data[i + 1][j]
      const west = data[i][j - 1]
      const east = data[i][j + 1]
      const northwest = data[i - 1][j - 1]
      const northeast = data[i - 1][j + 1]
      const southwest = data[i + 1][j - 1]
      const southeast = data[i + 1][j + 1]

      result[i][j] = (
        4 * center +
        2 * (north + south + west + east) +
        northwest + northeast + southwest + southeast
      ) / 16
    }
  }

  return result
}

平滑会减少噪声,也会改变局部峰谷值。展示型等高线可以接受,测绘或水文计算应保留原始网格,并明确平滑参数。

生成需要计算的高程层

根据最小高程、最大高程和等高距得到阈值:

const startHeight = Math.ceil(minHeight / interval) * interval
const endHeight = Math.floor(maxHeight / interval) * interval

for (
  let height = startHeight;
  height <= endHeight;
  height += interval
) {
  const segments = marchingSquares(
    data,
    height,
    expandedBounds,
  )
}

若网格为 rows × columns,共有 levelCount 个高程层,主要计算量近似为:

O(levelCount × rows × columns)

等高距越小,高程层越多。采样密度和等高距会同时放大计算量。

Marching Squares 如何提取线段

算法逐个检查网格单元。一个单元有四个角,把“高于阈值”记为 1,“低于阈值”记为 0,就能得到 0 到 15 的状态编号:

let caseIndex = 0
if (v0 >= threshold) caseIndex |= 1
if (v1 >= threshold) caseIndex |= 2
if (v2 >= threshold) caseIndex |= 4
if (v3 >= threshold) caseIndex |= 8

状态 0 表示四角都低于阈值,状态 15 表示四角都高于阈值。这两种情况没有等高线穿过单元。

其余状态需要计算等高线和四条边的交点。交点不是简单中点,而是根据两端高程线性插值:

function interpolateEdge(
  start: number,
  end: number,
  startValue: number,
  endValue: number,
  threshold: number,
) {
  if (startValue === endValue) return (start + end) / 2

  const ratio =
    (threshold - startValue) / (endValue - startValue)

  return Cesium.Math.lerp(start, end, ratio)
}

每个状态对应一条或两条线段:

switch (caseIndex) {
  case 1:
  case 14:
    return [[top, left]]
  case 2:
  case 13:
    return [[top, right]]
  case 3:
  case 12:
    return [[left, right]]
  case 4:
  case 11:
    return [[right, bottom]]
  case 6:
  case 9:
    return [[top, bottom]]
  case 7:
  case 8:
    return [[bottom, left]]
}

状态 5 和 10 是鞍点,同一个角点状态存在两种连接方式。可以用单元中心平均高程决定连接方向:

const centerAverage = (v0 + v1 + v2 + v3) / 4

if (caseIndex === 5) {
  return centerAverage >= threshold
    ? [[top, left], [bottom, right]]
    : [[top, right], [bottom, left]]
}

把零散线段连接成折线

Marching Squares 返回的是每个单元内的短线段。渲染前要按端点连接:

const keyOf = (point: LngLatPoint) =>
  `${point.x.toFixed(9)},${point.y.toFixed(9)}`

const endpointMap = new Map<string, number[]>()

segments.forEach((segment, index) => {
  for (const point of segment) {
    const key = keyOf(point)
    const list = endpointMap.get(key) ?? []
    list.push(index)
    endpointMap.set(key, list)
  }
})

从一条未使用线段出发,分别从头部和尾部查找相邻线段,直到找不到共享端点。toFixed(9) 是对浮点端点做量化,避免理论相同的坐标因为计算误差无法匹配。

裁剪回用户多边形

外接矩形中会生成多边形外部的等高线。每条连接后的折线还要执行:

  1. 用 AABB 快速排除完全不相交的线。
  2. 判断每个线段端点是否位于多边形内。
  3. 计算线段与多边形边界的交点。
  4. 按进出状态拆成若干内部折线。

简单的点在多边形内判断可以使用射线法:

function pointInPolygon(point: LngLatPoint, polygon: LngLatPoint[]) {
  let inside = false

  for (let i = 0, j = polygon.length - 1; i < polygon.length; j = i++) {
    const a = polygon[i]
    const b = polygon[j]
    const intersects =
      a.y > point.y !== b.y > point.y &&
      point.x < (b.x - a.x) * (point.y - a.y) / (b.y - a.y) + a.x

    if (intersects) inside = !inside
  }

  return inside
}

批量渲染而不是一条线一个 Entity

贴地等高线可以使用 GroundPolylineGeometry,把多个 GeometryInstance 合并到同一个 GroundPolylinePrimitive:

const instances = lines.map(line =>
  new Cesium.GeometryInstance({
    geometry: new Cesium.GroundPolylineGeometry({
      positions: line.map(point =>
        Cesium.Cartesian3.fromDegrees(point.x, point.y),
      ),
      width: 2,
    }),
    attributes: {
      color: Cesium.ColorGeometryInstanceAttribute.fromColor(
        Cesium.Color.CYAN,
      ),
    },
  }),
)

viewer.scene.primitives.add(
  new Cesium.GroundPolylinePrimitive({
    geometryInstances: instances,
    appearance: new Cesium.PolylineColorAppearance({ translucent: true }),
  }),
)

大量标签也应使用 LabelCollection,并限制数量。否则密集等高线会同时带来渲染负担和视觉噪声。

长循环可以每处理若干高程层后让出主线程:

if (layerIndex % 8 === 0) {
  await new Promise(resolve => setTimeout(resolve, 0))
}

这不会减少总计算量,但能让浏览器处理输入和绘制进度状态。

用生成版本淘汰旧任务

用户修改等高距或重新绘制区域时,上一轮采样和计算可能仍在执行:

let generationToken = 0

async function generate() {
  const token = ++generationToken
  const isStale = () => token !== generationToken

  const grid = await sampleTerrainGrid()
  if (isStale()) return

  const contours = await generateContours(grid, isStale)
  if (isStale()) return

  renderContours(contours)
}

function clear() {
  generationToken += 1
  removeContourPrimitives()
}

这个判断保证旧结果不会覆盖当前区域。它不能取消已经发出的地形请求,若需要真正取消,还要在数据请求层增加可取消机制。

适用边界

当前方法适合局部地形展示,还存在这些边界:

  • 规则网格在经纬度空间中生成,高纬度或大范围区域会产生明显尺度畸变。
  • 一次高斯平滑会改变原始峰谷高程。
  • 跨越 ±180° 经线的区域需要先展开经度。
  • 极小等高距会增加高程层数量并放大噪声。
  • 结果精度不会超过输入地形服务本身的分辨率。

局部等高线的核心不是画线,而是把异步地形采样结果转换成稳定的规则网格,再通过 Marching Squares 提取、连接、裁剪和批量渲染。每一步都要保留自己的输入和输出,参数变化时才知道哪些阶段可以复用。

参考资料

相关文章

Cesium 区域等高线:地形采样与 Marching Squares

从交互绘制区域开始,完成规则网格采样、等高线提取、平滑、标注和高程统计

Cesium 三维风场:用 WebGL 2 GPU 粒子描绘全球气流

拆解案例中的风场数据、GPU 粒子更新、轨迹渲染、颜色映射和 Vue 参数控制