【WEBGL】谷歌瓦片图加载从原理到实现

年前为 NothingJs 实现了一个扩展 NJ_lod_ground，目标是简单实现加载谷歌瓦片。为了让读者更加容易的理解，我直接改成了 WebGL 实例（总代码800行左右、依赖glMatrix）。并且把相关内容整理到本文，希望能帮助到刚刚入门的同学。

工程地址在文章结尾。

WGS84 大地坐标系 和 Web 墨卡托投影

GIS 领域最离不开的就是坐标变换，首先要搞清楚的就是地球上的一个点如何变换成地图上的一个点。文章不会详细讲解变换方法，因为本文重点并不是算法。但是还是要说清楚整个过程，我们已经清楚地球本身不是一个规则的球体，为了计算方便，需要有一个标准的大地坐标系来简化计算，而 WGS84（World Geodetic System一1984 Coordinate System） 就是这样一个坐标系。但是大地坐标系是三维坐标系，要映射到二维地图上还需要一步，就是投影变换（仿射变换），比如墨卡托投影。

Web 墨卡托 定义的大地坐标系是 WGS84 坐标系，投影方式与墨卡托投影类似，但是投影时地球不再当做是椭球体而是半径是6378137米的标准球体。

首先我们先简单了解一下 Web 墨卡托投影 的历史：

• 2005年 - 谷歌在谷歌地图中首次使用，当时的 Web墨卡托 使用者还称其为 世界墨卡托（World Mercator）- Spherical Mercator (unofficial deprecated ESRI)，代号 WKID 54004。
• 2006年 - OSGeo 在提出的 TMS - Tile Map Service 标准中使用代号 OSGEO:41001（WGS84 / Simple Mercator）- Spherical Mercator (unofficial deprecated OSGEO / Tile Map Service)。
• 2007年 - Christopher Schmidt（OpenLayers的重要贡献者之一） 在通过一次 GIS 讨论中为了在 OpenLayers 中使用谷歌投影，提出给 Web墨卡托 使用一个统一的代号 900913 - 形似 Google，并在OpenLayers的 OpenLayers/Layer/SphericalMercator.js 中正式使用代号 900913。
• 2008年 - EPSG 在6.15版本中正式给谷歌地图投影赋予 CRS 代号 EPSG:3785（Popular Visualisation CRS / Mercator），这也是 Web墨卡托 正式被 EPSG 组织承认。
• 2009年 - EPSG 使用新代号 EPSG:3857 代替之前的 EPSG:3785，给谷歌地图投影方法命名为 “公共可视化伪墨卡托投影（PVPM）”。
• 至今 - EPSG:3857（WGS 84 / Pseudo-Mercator） 代号是 web墨卡托 的正式代号。

谷歌瓦片

经过投影变换后，地理坐标就变成了平面地图坐标。考虑到需要地图的精度有大有小（缩放），所以将地图分级：顶层为0级，由一张256像素见方的图片存储，向后每多一级，像素是当前级别的4倍。由此便组成了一个金字塔式的地图瓦片层级结构，每张瓦片的大小固定为256像素的方形。

有了地图瓦片还需要对地图瓦片进行编号才行，谷歌采用XYZ表示瓦片的坐标和显示级别（缩放级别），其中XY的原点在左上角，X从左向右，Y从上向下，Z则表示显示级别（缩放级别）。

假如我们需要一张256像素的世界卫星地图，我们可以在浏览器访问：http://mt2.google.cn/vt/lyrs=s&hl=zh-CN&gl=cn&x=0&y=0&z=0

如何决定缩放级别

了解了投影和瓦片获取方式，还是不能实现一个简单的地图，我们还需要知道什么条件下加载哪些瓦片图。我程序里处理方式比较简单，需要知道当前显示级别，然后根据显示级别和摄像机与地平面的交点推导出中心瓦片坐标，然后在按中心瓦片坐标计算当前需要加载的全部瓦片坐标。

最重要的一点就是，如何获取当前级别呢？我并没有去想计算级别的方法，因为 WebGIS 领域里，开源的 Cesium 发展非常不错，于是就在 Cesium 源码中搜索了一下，找到了相关代码，借鉴这部分代码很快就完成了我们需要的方法。

不多说，上代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

// 通过给定参数计算纹素间距, 核心就是理解视坐标
// 注意：这里为了计算方便 仅仅算了纵向的像素间距
function getCurrentTexelSpacing(camera, screen, pos) {
  // 相机 和 相机朝向与地面的交点 的距离
  const distanceToPos = glMatrix.vec3.distance(camera.position, pos);

  // 实际上是 
  // Math.tan(相机yFov / 2.0) * distanceToPos = 纹素间距 * 屏幕高度 / 2.0
  return (2.0 * Math.tan(camera.yfov / 2.0) * distanceToPos) / screen[1];
}

// 通过零级别的纹素间距定义当前级别
// 零级别的纹素间距是一定的，即 2.0 * Math.PI * 地球半径 / 256 = 地球周长 / 256
function getLevel(texelSpacing) {
  const twoToTheLevelPower = levelZeroTexelSpacing / texelSpacing;
  // 实际上计算的是 2 的几次方是 twoToTheLevelPower 这样就不难理解了
  const level = Math.log(twoToTheLevelPower) / Math.log(2);
  const rounded = Math.round(level);

  return Math.max(rounded | 0, 0);
}

// camera 相机
// screen 渲染的长宽
// pos 相机朝向与地面夹角（这个夹角算法也是 Cesium 里的原始方法）
function getCurrentLevel(camera, screen, pos) {
  const texelSpacing = getCurrentTexelSpacing(camera, screen, pos);

  return getLevel(texelSpacing);
}

核心代码就这些，getCurrentLevel 就是获取当前显示级别的方法（当然我对原始方法稍微修改了一下以适应本地代码）。

通过代码以及注释不难理解，实际上缩放级别是通过 显示级别0下的纹素间距 和 摄像机位置与其朝向同地面的交点间距离上的纹素间距 比值决定的。

坐标转换

知道了当前级别还是不能实现一个简单的地图，因为是使用 WebGL 绘制的，所以还需要知道如何通过级别确定需要绘制的顶点坐标，这个就涉及到坐标转换的问题了。

假如我们知道了一个瓦片的 XYZ，如何获取一个瓦片所涉及到的正方形在世界坐标系下的顶点信息呢？有了目标我们再查看 Cesium ，就能很快找到需要的代码，然后借鉴之。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

// 获取某个级别横向瓦片数量
function getNumberOfXTilesAtLevel(level) {
  return numberOfLevelZeroTilesX << level;
}

// 获取某个级别纵向瓦片数量
function getNumberOfYTilesAtLevel(level) {
  return numberOfLevelZeroTilesY << level;
}

// 瓦片坐标 转 经纬度边界坐标（东西南北）
function tileXYToRadians(x, y, level) {
  // 获取整个世界地图的区域边界（弧度）
  const rectangle = boundary;
  // 获取当前显示级别横向纵向瓦片数量
  const xTiles = getNumberOfXTilesAtLevel(level);
  const yTiles = getNumberOfYTilesAtLevel(level);

  // 计算每个瓦片的宽度（弧度）
  const xTileWidth = rectangle.width / xTiles;
  const west = x * xTileWidth + rectangle.west;
  const east = (x + 1) * xTileWidth + rectangle.west;

  // 计算每个瓦片的高度（弧度）
  const yTileHeight = rectangle.height / yTiles;
  const north = rectangle.north - y * yTileHeight;
  const south = rectangle.north - (y + 1) * yTileHeight;

  // 返回边界信息 弧度
  return {
    west,
    south,
    east,
    north
  };
}

// 弧度经纬度信息 转换成 世界坐标
function tileLonlatRadians2world(lonlat) {
  const x = ((lonlat[0]) / boundary.width) * WIDTH;
  const y = ((-lonlat[1]) / boundary.height) * HEIGHT;

  return [x, y];
}

由源码可知，通过 tileXYToRadians 可以根据当前瓦片信息转换成边界弧度信息，再通过 tileLonlatRadians2world 由弧度信息可以获取世界坐标信息。

预备知识

好了，我们了解了投影，瓦片，级别和坐标之后，我们就可以着手写这个瓦片加载程序了。

但是再看源码之前还需要了解一些相关知识：

• WebGL API。
• 前端鼠标事件相关操作，可以查看源码了解。
• 前端类的封装以及闭包方式（为了简化开发环境和提高兼容性，源码并没有采用es6编写）。

工程文件说明

实例

已知问题

这部分是因为工程本身是为了简化代码说明瓦片加载过程，所以就忽略了部分问题，如果读者有兴趣可以自己着手解决这部分问题。

• 采样 - 目前使用的是前向渲染，用的是默认硬件采样方式。
• 坐标抖动 - 这个是坐标计算精度的问题。
• 摄像机操作 - 摄像机操作还比较生硬，平移速度没有根据缩放而变化；嵌入iframe鼠标操作并没有兼容。
• 瓦片加载方式 - 瓦片加载方式是按级别加载的， 并不会平滑过渡，也不会按照需要加载不同级别的瓦片。

引用

• GCJ02/BD09/80/54/2000常用坐标系详解-第六节（二）
• 瓦片地图原理
• WebGIS前端瓦片地图显示原理及实现
• 从底层谈，WebGIS 原理、设计、实现
• Web墨卡托投影
• Web Mercator投影与反投影计算公式
• 墨卡托投影坐标系（Mercator Projection）原理及实现C代码
• OpenStreetMap/Google/百度/Bing瓦片地图服务(TMS)

源码