How I built a wind map with WebGL

[XXHolic]

目录

• 引子
• 正文
• 基于 CPU 的风场可视化
• OpenGL 基础
• 获取风场数据
• 基于 GPU 移动粒子
• 绘制粒子
• 绘制粒子轨迹
• 风场插值查找
• GPU 上的伪随机生成器
• 下一步是什么？
• 参考资料

引子

对风场可视化的效果感兴趣，搜资料的时候发现这篇文章，读了后觉得翻译一下以便再次查阅。

• 原文：How I built a wind map with WebGL

正文

查看我的基于 WebGL 的风场模拟演示！让我们深入了解它的工作原理。

我要坦白的说：在 Mapbox 工作的最后几年里，我像躲避瘟疫一样避免直接的 OpenGL/WebGL 编程。原因之一：OpenGL API 和术语让我深感恐惧。它看起来总是那么复杂、混乱、丑陋和冗长，以至于我永远都无法投入其中。只要听到 stencil masks、 mipmap、depth culling、 blend functions、 normal maps 等术语，我就会有一种不安的感觉。

今年，我最终决定直面我的恐惧，使用 WebGL 构建一些有意义的东西。2D 风场模拟看起来是一个完美的机会——它很有用，视觉上令人惊叹，而且具有挑战性，但在能力范围内感觉仍然可以实现。我惊讶地发现，它远没有看上去那么可怕！

基于 CPU 的风场可视化

网上有很多风场可视化的例子，但最受欢迎和最有影响的是 Cameron Beccario 的著名项目 earth.nullschool.net 。它本身不是开源的，但它有一个旧的开源版本，大多数其它实现都是基于这个版本编写代码的。一个著名的开源派生是 Esri Wind JS 。使用该技术的流行气象服务包括 Windy 和 VentuSky 。

通常，浏览器中的这种可视化依赖于 Canvas 2D API，大致如下所示：

1. 在屏幕上生成一组随机粒子位置并绘制粒子。
2. 对于每个粒子，查询风的数据以获取其当前位置的粒子速度，并相应地移动它。
3. 将一小部分粒子重置到随机位置。这样可以确保风吹走的区域永远不会完全变空。
4. 逐渐淡出当前屏幕，并在顶部绘制新定位的粒子。

这样做会有随之而来的性能限制：

• 风粒子的数量需保持较低（例如，地球示例使用~5k）。
• 每次更新数据或视图时都会有很大的延迟（例如，地球示例大约 2 秒），因为数据处理成本很高，而且发生在 CPU 端。

此外，要将其集成为基于 WebGL 的交互式地图（如 Mapbox）的一部分，你必须在每一帧上将画布元素的像素内容加载到 GPU ，这将大大降低性能。

我一直在寻找一种方法，用 WebGL 在 GPU 端重新实现完整的逻辑，这样它会很快，能够绘制数百万个粒子，并且可以集成到 Mapbox GL 地图中，而不会造成很大的性能损失。幸运的是，我偶然发现了 Chris Wellons 所写关于 WebGL 中粒子物理 的精彩教程，并意识到风场可视化可以使用相同的方法。

OpenGL 基础

令人困惑的 API 和术语使得 OpenGL 图形编程非常难以学习，但从表面上看，这个概念非常简单。这里有一个实用的定义：

OpenGL 为高效绘制三角形提供了 2D API。

所以基本上你用 GL 所做的就是画三角形。除了可怕的 API 之外，困难还来自于执行此操作所需的各种数学和算法。它还可以绘制点和基本线（无平滑或圆形连接/封口），但很少使用。

OpenGL 提供了一种特殊的类 C 语言—— GLSL ——来编写由 GPU 直接执行的程序。每个程序分为两部分，称为着色器——顶点着色器和片元着色器。

顶点着色器提供用于转换坐标的代码。例如，将三角形坐标乘以 2 ，使我们的三角形看起来两倍大。我们在绘图时传递给 OpenGL 的每个坐标都将运行一次。一个基本的例子：

attribute vec2 coord;
void main() {
    gl_Position = vec4(2.0 * coord, 0, 1);
}

片元着色器提供用于确定每个绘制像素颜色的代码。你可以用它做很多很酷的数学运算，但最终它类似“把三角形的当前像素画成绿色”。示例：

void main() {
    gl_FragColor = vec4(0, 1, 0, 1);
}

在顶点着色器和片元着色器中，都可以做的一件很酷的事情是添加一个图像（称为纹理）作为参数，然后在该图像的任何点中查找像素颜色。我们将在风场可视化中很依赖这个。

片元着色器代码的执行是大规模并行的，并且硬件加速很快，因此通常比 CPU 上的等效计算快很多数量级。

获取风场数据

美国国家气象局每 6 小时发布一次全球天气数据，称为 GFS，以纬度/经度网格的形式发布相关数值（包括风速）。它以一种称为 GRIB 的特殊二进制格式编码，可以使用一组特殊的工具将其解析为人类可读的 JSON 。

我编写了几个小脚本，下载并将风数据转换成一个简单的 PNG 图像，风速编码为 RGB 颜色——每个像素的水平速度为红色，垂直速度为绿色。看起来是这样的：

你可以下载更高分辨率的版本（2x和4x），但 360×180 网格对于低缩放可视化来说已经足够了。PNG 压缩非常适合这种数据，上面的图像通常只有 80 KB 左右。

基于 GPU 移动粒子

现有风场可视化将粒子状态存储在 JavaScript 数组中。我们如何在 GPU 端存储和操作该状态？一种称为计算着色器（在 OpenGL ES 3.1 和等效的 WebGL 2.0 规范中）的新 GL 功能允许你在任意数据上运行着色器代码（无需任何渲染）。不幸的是，跨浏览器和移动设备对新规范的支持非常有限，因此我们只剩下一个实用选项：纹理。

OpenGL 不仅允许你绘制屏幕，还允许绘制纹理（通过称为帧缓冲区的概念）。因此，我们可以将粒子位置编码为图像的 RGBA 颜色，将其加载到 GPU ，在片着色器中根据风速计算新位置，将其重新编码为 RGBA 颜色，并将其绘制到新图像中。

X 和 Y 为了存储足够的精度，我们将每个组件存储在两个字节中——分别为 RG 和 BA ，为每个组件提供 65536 个不同值的范围。

一个 500×500 的示例图像将容纳 250000 个粒子，我们将使用片元着色器移动每个粒子。生成的图像如下所示：

以下是在片元着色器中如何从 RGBA 中解码和编码位置：

// lookup particle pixel color
vec4 color = texture2D(u_particles, v_tex_pos);
// decode particle position (x, y) from pixel RGBA color
vec2 pos = vec2(
    color.r / 255.0 + color.b,
    color.g / 255.0 + color.a);
... // move the position
// encode the position back into RGBA
gl_FragColor = vec4(
    fract(pos * 255.0),
    floor(pos * 255.0) / 255.0);

在下一帧中，我们可以将这个新图像作为当前状态，并将新状态绘制到另一个图像中，以此类推，每帧交换两个状态。因此，借助两个粒子状态纹理，我们可以将所有风模拟逻辑移动到 GPU 。

这种方法的速度非常快，我们不需要在浏览器上每秒更新 5000 个粒子 60 次，而是可以突然处理一百万个。

需要记住的一点是，在两极附近，粒子沿 X 轴的移动速度应该比赤道上的粒子快得多，因为相同的经度表示的距离要小得多。以下着色器代码可以处理这一点：

float distortion = cos(radians(pos.y * 180.0 - 90.0));
// move the particle by (velocity.x / distortion, velocity.y)

绘制粒子

正如我前面提到的，除了三角形，我们还可以绘制基本的点——很少使用，但非常适合这样的 1 像素粒子。

要绘制每个粒子，我们只需在顶点着色器中的粒子状态纹理上查找其像素颜色以确定其位置；然后通过从风纹理查找其当前速度来确定片元着色器中的粒子颜色；最后将其映射到一个漂亮的颜色渐变（我从可靠的 ColorBrewer2 中选择颜色）。在这一点上，它看起来是这样的：

如果有点空隙，那是一些东西。但单凭粒子运动很难获得风向感。我们需要添加粒子轨迹。

绘制粒子轨迹

我尝试的第一种绘制轨迹的方法是使用 WebGL 的 PreserveDrawingBuffer 选项，它使屏幕状态在帧之间保持不变，这样我们可以在粒子移动时在每一帧上反复绘制粒子。然而，这个 WebGL 特性是一个巨大的性能打击，许多文章建议不要使用它。

相反，与使用粒子状态纹理的方式类似，我们可以将粒子绘制到纹理中（该纹理依次绘制到屏幕上），然后在下一帧上使用该纹理作为背景（稍微变暗），并每一帧交换输入/目标纹理。除了更好的性能之外，这种方法的一个优点是我们可以将其直接移植到本机代码（没有与 preserveDrawingBuffer 等效的代码）。

风场插值查找

在纬度/经度栅格上，风数据针对特定点有对应的值，例如（50,30）、（51,30）、（50,31）、（51,31）地理点。如何获得任意中间值，例如（50.123,30.744）？

OpenGL 在查找纹理颜色时提供自带插值。然而，它仍然会导致块状、像素化的图案。以下是在缩放时，在风纹理中这些瑕疵的示例：

幸运的是，我们可以通过在每个风探测器中查找 4 个相邻像素，并在片元着色器中的本地像素上对其进行手动双线性插值计算，来平滑瑕疵。它的成本更高，但修复了瑕疵并产生更流畅的风场可视化。以下是与此技术相同的区域：

GPU 上的伪随机生成器

还有一个棘手的逻辑需要在 GPU 上实现——随机重置粒子位置。如果不这样做，即使是大量的风粒子也会变为屏幕上的几行，因为风吹走的区域会随着时间变空：

问题是着色器没有随机数生成器。我们如何随机决定粒子是否需要重置？

我在 StackOverflow 上找到了一个解决方案——一个用于生成伪随机数的 GLSL 函数，它接受一对数字作为输入：

float rand(const vec2 co) {
    float t = dot(vec2(12.9898, 78.233), co);
    return fract(sin(t) * (4375.85453 + t));
}

这个奇特的函数依赖于 sin 的结果变化。然后我们可以这样做：

if (rand(some_numbers) > 0.99)
    reset_particle_position();

这里的挑战在于为每个粒子选择一个足够“随机”的输入，以便生成的值在整个屏幕上是一致的，并且不会显示奇怪的图案。

使用当前粒子位置作为源并不完美，因为相同的粒子位置将始终生成相同的随机数，因此某些粒子将在同一区域消失。

使用在状态纹理中的粒子位置也不起作用，因为相同的粒子将始终消失。

我最终得到的结果取决于粒子位置和状态位置，再加上在每一帧上计算并传递给着色器的随机值：

vec2 seed = (pos + v_tex_pos) * u_rand_seed;

但我们还有另一个小问题——粒子速度非常快的区域看起来比没有太多风的区域密度要大得多。我们可以通过对更快的粒子增加粒子重置速率来平衡这一点：

float dropRate = u_drop_rate + speed_t * u_drop_rate_bump;

这里的 speed_t 是一个相对速度值（从0到1），u_drop_rate 和 u_drop_rate_bump 是可以在最终可视化中调整的参数。以下是它如何影响结果的示例：

下一步是什么？

结果是一个完全由 GPU 驱动的风场可视化，可以以 60fps 的速度渲染一百万个粒子。试着在演示中使用滑块，并查看最终的代码——总共大约 250 行，我努力使其尽可能的可读。

下一步是将其集成到可以探索的实时地图中。我在这方面取得了一些进展，但还不足以分享一个实时演示。这里有一些部分片段：

感谢你的阅读，请继续关注更多更新！如果你错过了，请查看我上一篇关于空间算法的文章。

非常感谢我的 Mapbox 队友 kkaefer 和 ansis，他们耐心地回答了我所有关于图形编程的傻问题，给了我很多宝贵的提示，帮助我学到了很多东西。❤️

Back to top ⬆️

参考资料

• How I built a wind map with WebGL

🗑️

最近听音乐的时候，发现了央视的节目《国家宝藏》，国乐演奏真是震撼！

• 《水龙吟》《九州同》《醉太平》《象王行》《破阵乐》合集