GeoTIFF 高程数据可视化：从不可见到可感知
前言
一张 GeoTIFF 高程图片，每个像素存储的不是颜色，而是一个浮点数——代表该点的海拔高度。这样的图片无法直接展示给人眼，甚至很多 TIFF 都只有单通道。那么，如何让这些「看不见」的数据变成「看得懂」的地形图呢？
本文将介绍两个核心技术：山体阴影（Hillshade） 和 色带映射（Color Ramp） ，并结合 WebGL2 GPU 加速，实现高程数据的实时可视化。
1. 问题：为什么 TIFF 高程图无法直接展示？
普通图片的每个像素是 RGB 颜色值（0~255），人眼可以直接辨识。而高程 TIFF 中的像素值可能是：
像素[0][0] = 342.7 （海拔342.7米）
像素[0][1] = 343.1
像素[1][0] = 341.9
...
这些浮点数没有颜色含义，且数据范围远超 0~255。即使强制显示为灰度图，你看到的只是一片近乎均匀的灰色——因为相邻像素的高程差异太小，人眼根本分辨不出。
我们需要两种手段来让高程信息「可见」：
- 山体阴影：通过模拟光照产生明暗对比，让地形起伏变得可感知
- 色带映射：将高程值映射到一组有意义的颜色上，让海拔高低一目了然
2. 山体阴影（Hillshade）的原理
山体阴影的核心思路是：如果知道地表每一点的朝向（法线），就能计算出阳光照射到该点时的明暗程度。
兰伯特余弦定理Lambert’s cosine law
一个表面接收的光照能量与 法线 n 和 光照方向 l 的余弦相关。见下图。
我们现在只有一张图片，里面的每个像素表示一个高程。那么如何获取地表的法线呢？
2.1 从高程到梯度：Sobel 卷积
要知道地表的朝向，首先需要知道地形在 x 和 y 方向上的变化率（梯度）。这通过 Sobel 算子 对高程数据做 3×3 邻域卷积来实现。
对于当前像素 C，采样其周围的 8 个像素：
tl tp tr
l [C] r
bl b br
Sobel 算子对 x 方向和 y 方向分别定义了权重矩阵：
Gx = [[-1, 0, 1], Gy = [[ 1, 2, 1],
[-2, 0, 2], [ 0, 0, 0],
[-1, 0, 1]] [-1, -2, -1]]
对一个像素周围的所有的像素取出，获取它们的高程值，在通过上述权重矩阵进行加权平均。
计算梯度代码如下：
dz/dx = ((tr + 2r + br) - (tl + 2l + bl)) / (8 × cellSizeX)
dz/dy = ((tl + 2tp + tr) - (bl + 2b + br)) / (8 × cellSizeY)
dz/dx 表示高程对x轴的变化速度 dz/dy 表示高程对y轴的变化速度
其中
cellSizeX和cellSizeY是每个像素在实际地理空间中代表的距离（单位：米），这样梯度的物理意义就是真实坡度（米/米） 。
2.2 从梯度到法线
有了 x/y 方向的梯度后，地表法线向量可以直接构造：
normal = normalize(-dz/dx × zFactor, -dz/dy × zFactor, 1.0)
zFactor是高程夸张系数，值越大，山体起伏在视觉上越明显- 法线的 z 分量为 1.0，因为梯度已经是真实坡度（米/米），无需再做缩放
- 取负号是因为梯度方向指向高程增大的方向，而法线应该指向地表的「外侧」
2.3 Lambert 漫反射光照
有了法线和光线方向，就可以用 Lambert 余弦定理 计算照明强度：
shade = max(dot(normal, lightDir), 0.0) × intensity
lightDir是归一化后的光线方向向量dot(normal, lightDir)即两个向量的点积，等于cos(θ)，θ 为法线与光线的夹角- 当面朝光线时 cos(θ) ≈ 1，完全背光时 cos(θ) ≤ 0（取 0）
intensity控制光照强度
这是最经典的漫反射模型：法线越朝向光源越亮，越背离光源越暗。
最终输出为灰度图：亮处代表阳面，暗处代表阴面，人眼就能直观地感受到山体的立体起伏。
3. 色带映射（Color Ramp）
山体阴影解决了「分辨起伏」的问题，但无法区分海拔高低。色带映射通过将高程值映射到一组预定义颜色来解决这一问题。
3.1 高程归一化
首先将高程值归一化到 [0, 1] 区间：
t = (elevation - minElevation) / (maxElevation - minElevation)
3.2 色带控制点
定义一组控制点，每个控制点指定一个归一化位置和对应的 RGBA 颜色：
const TERRAIN_COLOR_RAMP: ColorStop[] = [ { offset: 0.0, color: [0, 97, 0, 128] }, // 深绿 — 低海拔
{ offset: 0.15, color: [16, 122, 0, 128] }, // 绿
{ offset: 0.3, color: [132, 173, 54, 128] }, // 黄绿
{ offset: 0.45, color: [202, 204, 68, 128] }, // 黄
{ offset: 0.6, color: [185, 152, 90, 128] }, // 棕黄
{ offset: 0.75, color: [148, 107, 62, 128] }, // 棕
{ offset: 0.9, color: [178, 178, 178, 128] }, // 灰
{ offset: 1.0, color: [255, 255, 255, 128] }, // 白 — 高海拔/雪线
];
这是经典的 hypsometric tints（等高着色）方案：低处绿色（植被），中部棕黄（裸岩），高处灰白（雪线）。
3.3 插值生成 1D 纹理
将控制点之间进行线性插值，生成一条固定宽度（如 256 像素）的颜色带：
// 对每个纹素位置 t ∈ [0, 1]
// 找到 t 所在的两个相邻控制点 lower 和 upper
// 在两者之间做线性插值
factor = (t - lower.offset) / (upper.offset - lower.offset)
color = lower.color + (upper.color - lower.color) × factor
这个颜色带作为 1D 查找纹理上传到 GPU，片元着色器中只需一次 texture() 采样即可查到对应颜色。
4. WebGL2 GPU 加速实现
上述算法如果在 CPU 上逐像素计算，对于大尺寸的 TIFF 会非常慢。通过 WebGL2，我们将整个计算过程交给 GPU 并行执行。
4.1 整体流水线
高程数据 (Float32Array)
│
▼
┌──────────────────────────────────────────┐
│ GPU 纹理上传 │
│ ┌──────────────┐ ┌──────────────────┐ │
│ │ 高程纹理 R32F │ │ 色带纹理 RGBA8 │ │
│ │ 纹理单元 0 │ │ 纹理单元 1 │ │
│ └──────┬───────┘ └────────┬─────────┘ │
│ │ │ │
│ ▼ ▼ │
│ ┌──────────────────────────────────┐ │
│ │ 片元着色器（每像素并行执行） │ │
│ │ │ │
│ │ 1. NoData 检测 → 透明输出 │ │
│ │ 2. 高程归一化 → t ∈ [0, 1] │ │
│ │ 3. 色带纹理采样 → RGBA 颜色 │ │
│ │ 4. Sobel 3×3 → 梯度 dz/dx, dy │ │
│ │ 5. 梯度 → 地表法线 │ │
│ │ 6. Lambert 光照 → shade 值 │ │
│ │ 7. 环境光混合 │ │
│ │ 8. 色带颜色 × 阴影亮度 → 输出 │ │
│ └──────────────┬───────────────────┘ │
│ │ │
│ ▼ │
│ ┌──────────────────────────────────┐ │
│ │ readPixels → Uint8Array → ImageData │
│ └──────────────────────────────────┘ │
└──────────────────────────────────────────┘
4.2 关键实现细节
高程纹理使用 R32F 格式
高程值是浮点数，使用 gl.R32F（单通道 32 位浮点）格式存储，采样方式设为 NEAREST，避免高程值被纹理插值篡改：
gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32F, width, height, 0,
gl.RED, gl.FLOAT, heightData);
色带纹理使用 LINEAR 采样
色带纹理使用 LINEAR 采样，这样即使纹理分辨率不高，颜色过渡也是平滑的：
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.LINEAR);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.LINEAR);
cellSize：像素到真实距离的换算
地理坐标系中，经度方向的实际距离随纬度变化。取影像中心纬度进行近似：
const metersPerDegLat = 111320;
const metersPerDegLng = 111320 * Math.cos(centerLat × π / 180);
const cellSizeX = (经度跨度 × metersPerDegLng) / 图像宽度;
const cellSizeY = (纬度跨度 × metersPerDegLat) / 图像高度;
这保证了 Sobel 梯度的物理单位正确，山体阴影不会因为像素比例失真。
环境光（Ambient Light）防止全黑
纯 Lambert 模型下，完全背光的区域会变成纯黑。加入环境光的最低亮度：
float finalShade = mix(u_ambientLight, 1.0, shade);
当 shade = 0（完全背光）时，输出为 ambientLight（如 0.15），而非 0。
4.3 片元着色器核心代码
const fsSource = /* glsl */ `#version 300 es
precision highp float;
uniform sampler2D u_heightMap; // R32F 高程纹理（纹理单元 0）
uniform sampler2D u_colorRamp; // RGBA8 色带查找纹理（纹理单元 1）
uniform vec2 u_resolution; // (width, height)
uniform vec3 u_lightDir; // 归一化后的光线方向
uniform float u_intensity; // 光照强度
uniform float u_zFactor; // 高程夸张系数
uniform float u_minElevation; // 有效高程最小值
uniform float u_maxElevation; // 有效高程最大值
uniform float u_noDataValue; // NoData 标记值
uniform float u_hasNoData; // 是否启用 NoData 判断（1.0=启用, 0.0=禁用）
uniform float u_ambientLight; // 环境光最低亮度
uniform vec2 u_cellSize; // 每个像素在 x/y 方向代表的实际距离（米）
in vec2 v_texCoord;
out vec4 fragColor;
float getHeight(vec2 uv) {
return texture(u_heightMap, uv).r;
}
void main() {
float elev = getHeight(v_texCoord);
// ── NoData 检测：无效像素输出全透明 ──
if (u_hasNoData > 0.5 && elev == u_noDataValue) {
fragColor = vec4(0.0, 0.0, 0.0, 0.0);
return;
}
// ── 高程归一化到 [0, 1]，用于色带查找 ──
float elevRange = u_maxElevation - u_minElevation;
float t = (elevRange > 0.0)
? clamp((elev - u_minElevation) / elevRange, 0.0, 1.0)
: 0.5;
// ── 色带查找：利用 1D 纹理采样获取颜色（含透明度） ──
vec4 rampColor = texture(u_colorRamp, vec2(t, 0.5));
// ── Hillshade 阴影计算（Sobel 3x3 + Lambert） ──
vec2 texelSize = 1.0 / u_resolution;
float tl = getHeight(v_texCoord + vec2(-texelSize.x, texelSize.y));
float tp = getHeight(v_texCoord + vec2( 0.0, texelSize.y));
float tr = getHeight(v_texCoord + vec2( texelSize.x, texelSize.y));
float l = getHeight(v_texCoord + vec2(-texelSize.x, 0.0));
float r = getHeight(v_texCoord + vec2( texelSize.x, 0.0));
float bl = getHeight(v_texCoord + vec2(-texelSize.x, -texelSize.y));
float b = getHeight(v_texCoord + vec2( 0.0, -texelSize.y));
float br = getHeight(v_texCoord + vec2( texelSize.x, -texelSize.y));
// Sobel 梯度除以像素实际地理距离，得到真实坡度（米/米）
float dzdx = ((tr + 2.0 * r + br) - (tl + 2.0 * l + bl)) / (8.0 * u_cellSize.x);
float dzdy = ((tl + 2.0 * tp + tr) - (bl + 2.0 * b + br)) / (8.0 * u_cellSize.y);
// 地表法线（梯度已是真实坡度，z 分量为 1.0 即可）
vec3 normal = normalize(vec3(-dzdx * u_zFactor, -dzdy * u_zFactor, 1.0));
// Lambert 漫反射
float shade = max(dot(normal, u_lightDir), 0.0) * u_intensity;
shade = clamp(shade, 0.0, 1.0);
// ── 混合环境光，防止阴影面全黑 ──
float finalShade = mix(u_ambientLight, 1.0, shade);
// ── 最终输出：色带颜色 × 阴影亮度，保留色带透明度 ──
fragColor = vec4(rampColor.rgb * finalShade, rampColor.a);
}
`;
一次 Pass，所有计算在片元着色器中同时完成——无需多趟渲染或中间缓冲区。
我们通过 webgl的readPixels 获取图片像素数据，再把它贴图到地图上，最终效果如下所示：
观察到地形的细节还是不错的。
5. 完整数据流总结
GeoTIFF 文件
│
▼
┌─────────────────────┐
│ geotiff.js 解码 │
│ → 宽、高、通道数 │
│ → Float32 栅格数据 │
└─────────┬───────────┘
│
┌─────────┴───────────┐
│ │
▼ ▼
┌────────────────┐ ┌────────────────┐
│ 元数据统计 │ │ 色带控制点 │
│ min/max 高程 │ │ 线性插值 │
│ NoData 处理 │ │ → 1D RGBA 纹理 │
└───────┬────────┘ └───────┬────────┘
│ │
▼ ▼
┌─────────────────────────────────────┐
│ WebGL2 单 Pass 渲染 │
│ │
│ 高程纹理(R32F) + 色带纹理(RGBA8) │
│ ↓ │
│ 片元着色器： │
│ NoData → Sobel → 法线 → Lambert │
│ → 色带采样 → 混合 → 输出 │
└──────────────┬──────────────────────┘
│
▼
┌──────────┐
│ ImageData │
│ RGBA 可视化│
└──────────┘
6. 效果调参指南
| 参数 | 作用 | 建议值 |
|---|---|---|
lightSource.direction | 光线方向向量 (x, y, z) | [-1, -1, 2] 模拟西北方向阳光 |
lightSource.intensity | 光照强度 | 0.8 ~ 1.5 |
zFactor | 高程夸张系数 | 平原地区 2 |
ambientLight | 环境光最低亮度 | 0.1 ~ 0.2 |
colorRamp | 自定义色带 | 按需定义，支持任意 RGBA 控制点 |
总结
GeoTIFF 高程可视化的本质是两个映射的叠加：
- 高程 → 颜色（色带映射）：让人一眼区分海拔高低
- 高程 → 明暗（山体阴影）：让人感知地形的立体起伏
两者结合——色带颜色乘以阴影亮度——就得到了既有色彩信息又有立体感的地形可视化图像。借助 WebGL2 的并行计算能力，整个过程在一次 GPU Pass 中完成，即使面对大尺寸高程数据也能高效处理。
