文章目录
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) 是对浮点端点做量化,避免理论相同的坐标因为计算误差无法匹配。
裁剪回用户多边形
外接矩形中会生成多边形外部的等高线。每条连接后的折线还要执行:
- 用 AABB 快速排除完全不相交的线。
- 判断每个线段端点是否位于多边形内。
- 计算线段与多边形边界的交点。
- 按进出状态拆成若干内部折线。
简单的点在多边形内判断可以使用射线法:
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 提取、连接、裁剪和批量渲染。每一步都要保留自己的输入和输出,参数变化时才知道哪些阶段可以复用。