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

原创
2020/02/22 21:34
阅读数 1.9W

年前为 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编写)。

工程文件说明

路径 说明
css/main.css 全局样式
js/Camera.js 摄像机对象
js/gl-matrix-min.js 数学计算库压缩版
js/gl-matrix.js 数学计算库
js/mouse-utils.js 鼠标事件工具类
js/Quad.js 瓦片对象
js/tiles-utils.js 瓦片相关工具类
js/webgl-utils.js gl相关工具类
index.html html页面
main.js 主要运行过程

实例

地址

已知问题

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

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

引用

源码

地址

展开阅读全文
加载中

作者的其它热门文章

打赏
0
0 收藏
分享
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部