WebGIS 瓦片地图引擎实现之——地图瓦片加载计算原理介绍

ThinkGIS GIS 背景在蚂蚁搞空间数据可视化，L7地理可视化引擎作者 65 人赞同了该文章

. 背景

1.1 地图瓦片之前

在地图瓦片技术使用之前，用户使用在线地图，一般都是客户端把将要显示的地理范围传送到服务端，服务器端将地理范围内的地理数据都查询出来，然后在服务端按照预先定义的专题地图的配置样式渲染成地图图片，再返回给客户端浏览，这也就是 OGC（地理信息联盟）标准的 WMS（Web Map Service）服务。

在 Web 时代 WMS 服务有几个大的缺点，由于用户浏览的地理范围是不确定性质的，获取范围内的数据是不确定，有时候数据会很大，数据从服务端的数据库到客户端的过程中，服务器 IO 操作和网络传输就很耗时；服务端获取数据后会进行数据渲染出图，占用大量 CPU 资源；用户频繁操作地图，服务端出图传输给客户端浏览过程中耗时非常长。

1.2 什么是地图瓦片

为了解决这些问题，谷歌地图率先提出了 TMS（Tile Map Service）服务，预先在服务端分层（金字塔模型）切片全量渲染，提出了 Web-Mercator 投影，按照地图不同的显示级别切成了瓦片坐标的图片。

用户每次访问时候，根据当前的地理范围映射到瓦片坐标的图片索引（XYZ），然后从后端请求这些图片索引，客户端拿到图片后按照顺序依次渲染图片即可，整个过程中用户体验有明显的提升。很快 TMS 服务成为了 WebGIS 工业标准，在该技术的推动下，OGC（地理信息联盟）也发布了基于 TMS 的 WMTS 服务规范，国内外很多地图厂商也基于该技术生产了自己的切片地图服务。

1.3 栅格瓦片与矢量瓦片

栅格瓦片就是上述所说的地图瓦片，只是瓦片的数据形式是图片，以矢量数据组织的瓦片是矢量瓦片。栅格瓦片技术使用很广泛，也存在一些问。瓦片预切图时间长，数据更新不方便，按金字塔模型进行切图很耗时，每进入下一级别图片成 2^zoom 指数级别增加；服务端资源要求高，切片过程需要极大服务端渲染和切图的计算资源，存储数据冗余，在存储了原始数据外，还存储了冗余的切片数据。栅格瓦片地图样式单一，因为最终渲染是图片数据，不能满足自定义地图风格。

因上面这些问题 Mapbox 提出了 MVT 矢量瓦片切片技术，地图定制化极为方便，该技术目前最受欢迎，除二维地理数据切片外三维数据也深受此技术的影响诞生了 3DTiles 。当然说到这里这里瓦片预切图时间长与数据存储冗余并没有完全得倒解决，这里不再延续展开目前有的方案（数据动态切片与不切片）。

2. 瓦片地图生产

为了更好的使用瓦片服务，我们需要了解瓦片地图的生产过程，这里我们以栅格瓦片为例。下面分别以重要节点来简述这几个步骤。

2.1 数据投影及配图

2.1.1 数据投影

数据库中地理数据存储的形式一般都是大地坐标系（以纬度、经度和高度表示坐标元素），一般使用经纬度描述地球上点的位置，但是对于地图地理数据在二维平面内展示的场景，需要通过投影的方式将三维空间中的点映射到二维空间中，地图投影需要建立地球表面点与投影平面点的一一映射关系，映射方法有很多种，各有各的优缺点，互连网地图中常用 Web 墨卡托投影，它是墨卡托投影的变体。

Web 墨卡托投影是一个球面墨卡托的投影，它具有与球面墨卡托相同的属性，子午线都是等距的垂直线，角度在局部是正确的，而面积会随着离赤道越来越远而膨胀，以至于极地区域被严重夸大。由于 Web 墨卡托投影指定了在 WGS 84 椭球面模型上给出的测量坐标，但在投影的时候定义在球面上的，因此会有偏离，但在小比例尺民用场景精度可忽略不计，它主要的优点还是计算简单，数学运算复杂度低。

在缩放层级为 0 时，投影的过程就是将球面投影到一个正方形的平面上，这个平面就是世界坐标调整为左上角为（0, 0），右下角为（256, 256）的正方形，假设单位为像数，那地图投影在一个 256px * 256px 的图片上。

因为将极点投影在无穷远处，所以使用 Web墨卡托投影的地图无法显示极点，数据覆盖范围在经度（-180°～180°），纬度（-85.051129°～85.051129°）之间，投影的地图也是正方形。

2.1.2 地图配图

对数据进行投影之后，就需要根据不同的比例尺级别对数据进行配置，不同比例尺级别，显示的内容详略是不一样的，对数据进行分层级进行配图，比如层级大于 18 的时候，才显示一些街道级别的道路和 POI 等。再根据数据的不同类型进行样式配置，比如海洋填充显示为为淡蓝色，铁路显示为黑白相间的线段等。

数据的分层配图比较简单，但要配出一张好看的图需要处理好很多细节，专业的桌面软件有 ArcGIS、QGIS 等，感兴趣可体验简化在线版 Maputnik 配置矢量瓦片样式，这里不再展开。

2.2 地图切片

这里以谷歌地图规范的 TMS 为例，各地图厂商切片方式略有不一样，但都是基于 TMS 大同小异。

2.2.1 瓦片坐标系统

与谷歌地图相同的有高德地图、OpenStreetMap，一般都采用这种瓦片坐标系，将 Web 墨卡托投影地图的左上角作为瓦片坐标系起点，往左方向为 X 轴，X 轴与北纬 85.05º 重合且方向向左；往下方向为 Y 轴，Y 轴与东经 180º 重合且方向向下。瓦片坐标最小等级为 0 级，此时平面地图是一个像素为 256*256 的瓦片。在某一瓦片层级下，瓦片坐标的 X 轴和 Y 轴各有 2^zoom 个瓦片。

WMTS 规范中，Web 墨卡托投影地图的左下角作为瓦片坐标系起点，也就是坐标系原点为东经 180°，南纬 85.05°，X 轴向右，Y 轴向上，瓦片坐标最小等级也从 0 级开始，常见使用该规范的地图服务有天地图、腾讯地图、ArcGIS、超图等。

需要注意的是，微软的必应地图使用的编码规范，Z 的定义与谷歌相同，同一层级的瓦片不用 XY 两个维度表示，而只用一个整数表示，该整数服从四叉树编码规则；百度地图的，Z从 1 开始，在最高级就把地图分为四块瓦片，XY 的原点在经纬度为 0 的位置，X从左向右，Y从下向上。

2.2.2 坐标转换过程

在了解了数据投影后，我们需要将经纬度坐标转换为瓦片坐标，Web 墨卡托的投影计算公式如下图所示

公示中 λ 是用弧度表示的经度，而 φ 是用弧度表示的纬度。经纬度坐标转换为瓦片坐标其推导过程为经纬度 => 米 => 像素坐标 => 瓦片坐标

经纬度 => 米

将经纬度转为球面 Web 墨卡托投影的米制单位，投影到正方形

const LatLongToMeterXY = (longitude: number, latitude: number) => {
  // 地球的周长的一半 20037508.342789244 单位米
  const circumferenceHalf = Math.PI * 2 * 6378137 / 2.0
  const mx = longitude * circumferenceHalf / 180
  const temp = Math.log(Math.tan((90 + latitude) * (Math.PI / 360.0))) / (Math.PI / 180.0)
  const my = circumferenceHalf * temp / 180
  
  return { mx, my }
}

// [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]

米 => 像素坐标

const MetersToPixelXY = (mx: number, my: number, zoom: number, tileSize = 256) => {
  // 地球的周长的一半 20037508.342789244 单位米
  const circumferenceHalf = Math.PI * 2 * 6378137 / 2.0
  
  // 米/每像素
  const resolution = Math.PI * 2 * 6378137 / (tileSize * Math.pow(2, zoom))
  
  
  const px = (mx + circumferenceHalf) / resolution
  const py = (my + circumferenceHalf) / resolution
  
  return { px, py }
}

像素坐标 => 瓦片坐标

const PixelXYToTileXY = (px: number, py: number, tileSize = 256) => {
  const tx = Math.floor(px / tileSize);
  const ty = Math.floor(py / tileSize);
  
  return { tx, ty }
}

2.2.3 经纬度与瓦片坐标互转

根据上面的推导过程，很容易算出经纬度坐标转瓦片坐标，瓦片坐标转经纬度坐标的公式，转换公式如下：

经纬度 => 瓦片坐标

const LatLongToTileXY = (longitude: number, latitude: number, zoom: number) => {
  // 地球的周长的一半 20037508.342789244 单位米
  const circumferenceHalf = Math.PI * 6378137;
  const mx = (longitude * circumferenceHalf) / 180;
  const temp = Math.log(Math.tan((90 + latitude) * (Math.PI / 360.0))) / (Math.PI / 180.0);
  const my = (circumferenceHalf * temp) / 180;
  // 米/每像素
  const resolution = (Math.PI * 2 * 6378137) / (tileSize * Math.pow(2, zoom));
  
  // px = (mx + circumferenceHalf) / resolution
  const px = ((longitude + 180) / 360) * Math.pow(2, zoom) * tileSize;
  // py = (my + circumferenceHalf) / resolution
  const py = ((circumferenceHalf * temp) / 180 + circumferenceHalf) / (Math.PI * 2 * 6378137) / (tileSize * Math.pow(2, zoom));
  
  // tx = Math.floor(px / tileSize)
  const tx = Math.floor(((longitude + 180) / 360) * Math.pow(2, zoom));
  // ty = (1 - arsinh(tan(latitude * π / 180)) / π) * Math.pow(2, zoom)
  // 双曲正弦函数 arsinh(x) => ln(x + (x^2 + 1)^0.5)
  const ty = Math.floor(
    ((1 -
      Math.log(
        Math.tan((latitude * Math.PI) / 180) +
          1 / Math.cos((latitude * Math.PI) / 180)
      ) /
        Math.PI) /
      2) *
      Math.pow(2, zoom)
  );

  return { tx, ty };
};

瓦片坐标 => 经纬度

const TileXYToLatLong = (tx: number, ty: number, zoom: number) => {
  const longitude = (x / Math.pow(2, zoom)) * 360 - 180;
  
  const tmp = Math.PI - (2 * Math.PI * y) / Math.pow(2, zoom);
  
  // latitude = (180 / Math.PI) * Math.atan(sinh(tmp))
  // 双曲函数sinh(x)=(e^x - e^(-x)) * 0.5
  const latitude =  (180 / Math.PI) * Math.atan(0.5 * (Math.exp(tmp) - Math.exp(-tmp)));

  return { longitude, latitude };
};

2.3 瓦片地图服务发布

切片完成后的瓦片文件在一个大目录下，一般每个缩放层级是一个目录，每列是一个子目录，并且该列中的每个瓦片是一个文件，每个瓦片文件的路径名一般为缩放层级 + 列号 + 行号，比如 /z/x/y.png 这种形式，客户端通过拼接服务的域名加文件名来请求对应的瓦片。一般服务器考虑到并发请求原因会对地图瓦片的服务器会提供几个子域名，比如 OSM 地图服务：