绘制地图

之前一直很感兴趣在 WebGL 中如何展示地理信息，最近阅读了 deck.gl 的源码，学习到了很多相关知识。本文将介绍如何根据数据展示一个多边形，后续文章将介绍路径以及选中效果的实现。

数据格式

绘制基础图形不是难事，毕竟几何属性特征相对简单且固定。而对于一个不规则的多边形，我们首先需要描述它的所有顶点坐标。 GeoJSON 以 JSON 格式展现简单的地理信息。比如我们想描述一个多边形，coordinates 定义了全部顶点坐标（经纬度）：

{
    type: "FeatureCollection",
    features: [
    {
        type: "Feature",
        geometry: {
            type: "Polygon",
            coordinates: [
                [102.0, 0.0],
                [103.0, 1.0],
                [104.0, 0.0],
                [105.0, 1.0]
            ]
        }
    }
    ]
}

当然，除了多边形，地图中常见的路径，散点，甚至是更灵活的附加信息也都有对应的格式描述，具体可以参考 GeoJSON 的 TypeScript 定义文件，例如支持的全部几何类型：

// @types/geojson

/**
 * The valid values for the "type" property of GeoJSON geometry objects.
 * https://tools.ietf.org/html/rfc7946#section-1.4
 */
export type GeoJsonGeometryTypes = "Point" | "LineString" | "MultiPoint" | "Polygon" | "MultiLineString" |
    "MultiPolygon" | "GeometryCollection";

接下来我们只关心多边形这一种几何类型，对于原始的顶点数据我们需要做一些处理。

Tessellation & Triangulation

「Tessellation」是将一个表面分割成多个多边形的过程，虽然分割方式多种多样，但由于底层 API 以三角形为基础进行绘制，分割成多个基本三角形的过程也叫做「Triangulation」：

ear clipping 算法

那么如何分割三角形呢？一种常见的方法叫做「ear clipping」。如果一个顶点相邻的两个顶点连线不与多边形的任何一边相交，那么这个顶点就构成了「ear」。我们以下图左侧为例，在这个判断标准下，v2 v4 和 v5 就是「ear」。然后移除找到的 ear 如 v4，继续判定其相邻的顶点 v5 和 v3，其中 v5 构成了 ear。以此类推最终所有的 ear 都被移除，多边形也最终被分割成多个三角形。可参考论文 Triangulation By Ear Clipping。

在具体实现方面，Mapbox 开源了 earcut，其中判定 ear 的代码如下：

// check whether a polygon node forms a valid ear with adjacent nodes
function isEar(ear) {
    var a = ear.prev,
        b = ear,
        c = ear.next;

    if (area(a, b, c) >= 0) return false; // reflex, can't be an ear

    // now make sure we don't have other points inside the potential ear
    var p = ear.next.next;

    while (p !== ear.prev) {
        if (pointInTriangle(a.x, a.y, b.x, b.y, c.x, c.y, p.x, p.y) &&
            area(p.prev, p, p.next) >= 0) return false;
        p = p.next;
    }

    return true;
}

上述代码涉及到如何判断一个点是否在三角形中，在光栅化实现中一种常见的办法是将判定点和直线位置关系的「edge function」作用于三角形的三边。感兴趣的话可以阅读 a Practical Implementation (The Rasterization Stage)：

最终我们使用 earcut 得到了分割之后的 indices 数组，供后续 drawElements 使用。有一点需要注意，如果 indices 使用了 Uint32Array 存储，则需要使用对应的 gl.UNSIGNED_INT 类型，WebGL1 需要开启 OES_element_index_uint 扩展：

import earcut from 'earcut';

let indices = earcut([10,0, 0,50, 60,60, 70,10]); // returns [1,0,3, 3,2,1]
indices = new Uint32Array(indices);

const indexBuffer = gl.createBuffer();
gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, indexBuffer);
gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, indices, gl.STATIC_DRAW);
// WebGL1 中使用 gl.UNSIGNED_INT 需要额外开启扩展
gl.drawElements(gl.TRIANGLES, indices.length, gl.UNSIGNED_INT, 0);

现在我们已经可以将多边形分割成若干三角形了，接下来我们需要关心原始的顶点坐标数据。

Web Mercator 坐标系转换

我们的多边形顶点坐标是用 Web Mercator 坐标表示的，我们需要将其投影到最终的裁剪坐标才能使用，当然中途会经过世界坐标系的转换。

首先根据 Web Mercator 转换到世界坐标系公式：

参考 deck.gl 的 vertex shader 实现：

uniform float u_ProjectScale; // zoom level
uniform vec3 u_PixelsPerMeter;
uniform mat4 u_ModelMatrix;

const float TILE_SIZE = 512.0;
const float PI = 3.1415926536;
const float WORLD_SCALE = TILE_SIZE / (PI * 2.0);

vec2 project_mercator_(vec2 lnglat) {
    float x = lnglat.x;
    return vec2(
        radians(x) + PI,
        PI - log(tan(PI * 0.25 + radians(lnglat.y) * 0.5))
    );
}

vec4 project_position(vec4 position) {
    return u_ModelMatrix * vec4(
        project_mercator_(position.xy) * WORLD_SCALE * u_ProjectScale,
        project_scale(position.z),
        position.w
    );
}

最后通过 VP 矩阵投影到最终的裁剪坐标系：

attribute vec4 a_Position;
attribute vec4 a_Color;
uniform mat4 u_ViewProjectionMatrix;
varying vec4 v_Color;

vec4 project_to_clipspace(vec4 position) { 
    return u_ViewProjectionMatrix * position;
}      

void main() {
    v_Color = a_Color;

    vec4 worldPosition = project_position(a_Position);
    gl_Position = project_to_clipspace(worldPosition);
}

现在我们只剩下最后一个问题，那就是 vertex shader 中使用的 u_ViewProjectionMatrix 如何计算得到。

View 矩阵

learnwebgl 上的摄像机运动一节很清晰地介绍了摄像机的移动方式，即绕 u v n 三轴旋转或沿三轴移动：

根据 deck.gl 的定义，使用到以下三个变量，根据这三个变量我们就能计算出所需要的 View 矩阵：

pitch 绕 u 轴旋转，比如取 0 就是朝正下方观察，对应 “tilt” 动作
bearing：绕 n 轴旋转一定角度，比如取 0 的时候就是上北下南，对应 “cant” 动作
altitude：沿 n 轴移动，摄像机拉远拉近，对应 “dolly” 动作

export function getViewMatrix({
  pitch,
  bearing,
  altitude,
}) {
  const vm = createMat4();
  // 首先沿 n 轴移动
  mat4_translate(vm, vm, [0, 0, -altitude]);
  // 绕 u 轴旋转
  mat4_rotateX(vm, vm, -pitch * DEGREES_TO_RADIANS);
  // 绕 n 轴旋转
  mat4_rotateZ(vm, vm, bearing * DEGREES_TO_RADIANS);
  return vm;
}

接下来需要根据起始经纬度设置摄像机位置，转换公式在上面 vertex shader 中其实已经实现过一遍了，这里就不再展示了。唯一要注意的是 Web Mercator 坐标系和世界坐标系的 Y 轴方向是相反的：

// 起点经纬度坐标 -> 世界坐标
this.center = this._getCenterInWorld({longitude, latitude});

// 翻转 Y 轴方向
this.viewMatrixUncentered = mat4.scale([], viewMatrix, [1, -1, 1]);

// 摄像机位置移动到起点
this.viewMatrix = new Matrix4()
    .multiplyRight(this.viewMatrixUncentered)
    .translate(new Vector3(this.center || ZERO_VECTOR).negate());

至此我们就得到了完整的 View 矩阵，接下来是 Projection 矩阵的实现。

Projection 矩阵

投影矩阵没啥特别的，例如选择了透视投影：

new Matrix4().perspective({fovy, aspect, near, far});

根据 altitude 和 pitch，我们可以计算出以上投影矩阵需要的参数：

export function getProjectionParameters({
  width,
  height,
  altitude = DEFAULT_ALTITUDE,
  pitch = 0,
  nearZMultiplier = 1,
  farZMultiplier = 1
}) {
  // Find the distance from the center point to the center top
  // in altitude units using law of sines.
  const pitchRadians = pitch * DEGREES_TO_RADIANS;
  const halfFov = Math.atan(0.5 / altitude);
  const topHalfSurfaceDistance =
    Math.sin(halfFov) * altitude / Math.sin(Math.PI / 2 - pitchRadians - halfFov);

  // Calculate z value of the farthest fragment that should be rendered.
  const farZ = Math.cos(Math.PI / 2 - pitchRadians) * topHalfSurfaceDistance + altitude;

  return {
    fov: 2 * Math.atan((height / 2) / altitude),
    aspect: width / height,
    near: nearZMultiplier,
    far: farZ * farZMultiplier
  };
}

展示效果

现在我们已经能展示多边形，例如下面的美国地图，数据来自 Natural Earth：

Demo 地址

总结 & 后续

本文大量参考了 deck.gl 中的实现（尤其是 Uber 开源的 viewport-mercator-project 工具库），感兴趣的同学也可以直接阅读文档和源码。

在后续的文章中，我们将绘制地图中其他常见的几何图形，例如路径和散点，以及多边形描边（地理边界）。在常见的交互方面，摄像机移动以及基于 color-based picking 的多边形的拾取也都会介绍，任何问题都欢迎在评论区交流。

参考资料

deck.gl
GeoJSON
OpenStreetMap Wiki
Triangulation By Ear Clipping
RTR 3rd - 12.Polygonal Techniques

CATALOG

绘制多边形