从零打造一个Web地图引擎

Posted 2022-01-20 街角小林

说到地图，大家一定很熟悉，平时应该都使用过百度地图、高德地图、腾讯地图等，如果涉及到地图相关的开发需求，也有很多选择，比如前面的几个地图都会提供一套​​js API​​​，此外也有一些开源地图框架可以使用，比如​​OpenLayers​​​、​​Leaflet​​等。

那么大家有没有想过这些地图是怎么渲染出来的呢，为什么根据一个经纬度就能显示对应的地图呢，不知道没关系，本文会带各位从零实现一个简单的地图引擎，来帮助大家了解​​GIS​​基础知识及​​Web​​地图的实现原理。

选个经纬度

首先我们去高德地图上选个经纬度，作为我们后期的地图中心点，打开​​高德坐标拾取​​工具，随便选择一个点：

笔者选择了杭州的雷峰塔，经纬度为：​​[120.148732,30.231006]​​。

瓦片url分析

地图瓦片我们使用高德的在线瓦片，地址如下：

https://webrd01-4.is.autonavi.com/appmaptile?x=x&y=y&z=z&lang=zh_cn&size=1&scale=1&style=8

目前各大地图厂商的瓦片服务遵循的规则是有不同的：

谷歌XYZ规范：谷歌地图、OpenStreetMap、高德地图、geoq、天地图，坐标原点在左上角

TMS规范：腾讯地图，坐标原点在左下角

WMTS规范：原点在左上角，瓦片不是正方形，而是矩形，这个应该是官方标准

百度地图比较特立独行，投影、分辨率、坐标系都跟其他厂商不一样，原点在经纬度都为0的位置，也就是中间，向右为X正方向，向上为Y正方向

谷歌和​​TMS​​的瓦片区别可以通过该地址可视化的查看：​​地图瓦片​​。 虽然规范不同，但原理基本是一致的，都是把地球投影成一个巨大的正方形世界平面图，然后按照四叉树进行分层切割，比如第一层，只有一张瓦片，显示整个世界的信息，所以基本只能看到洲和海的名称和边界线，第二层，切割成四张瓦片，显示信息稍微多了一点，以此类推，就像一个金字塔一样，底层分辨率最高，显示的细节最多，瓦片数也最多，顶层分辨率最低，显示的信息很少，瓦片数量相对也最少：

每一层的瓦片数量计算公式：

Math.pow(Math.pow(2, n), 2)// 行*列：2^n * 2^n

十八层就需要​​68719476736​​张瓦片，所以一套地图瓦片整体数量是非常庞大的。

瓦片切好以后，通过行列号和缩放层级来保存，所以可以看到瓦片地址中有三个变量：​​x​​、​​y​​、​​z​​

x：行号
y：列号
z：分辨率，一般为0-18

通过这三个变量就可以定位到一张瓦片，比如下面这个地址，行号为​​109280​​，列号为​​53979​​，缩放层级为​​17​​：

https://webrd01.is.autonavi.com/appmaptile?x=109280&y=53979&z=17&lang=zh_cn&size=1&scale=1&style=8

对应的瓦片为：

关于瓦片的更多信息可以阅读​​瓦片地图原理​​。

坐标系简介

高德地图使用的是​​GCJ-02坐标系​​，也称火星坐标系，由中国国家测绘局在02年发布，是在GPS坐标（​​WGS-84​​坐标系）基础上经加密后而来，也就是增加了非线性的偏移，让你摸不准真实位置，为了国家安全，国内地图服务商都需要使用​​GCJ-02坐标系​​。

​​WGS-84​​坐标系是国际通用的标准，​​EPSG​​编号为​​EPSG:4326​​，通常GPS设备获取到的原始经纬度和国外的地图厂商使用的都是​​WGS-84​​坐标系。

这两种坐标系都是地理坐标系，球面坐标，单位为​​度​​，这种坐标方便在地球上定位，但是不方便展示和进行面积距离计算，我们印象中的地图都是平面的，所以就有了另外一种平面坐标系，平面坐标系是通过投影的方式从地理坐标系中转换过来，所以也称为投影坐标系，通常单位为​​米​​，投影坐标系根据投影方式的不同存在多种，在​​Web​​开发的场景里通常使用的是​​Web墨卡托投影​​，编号为​​EPSG:3857​​，它基于​​墨卡托投影​​，把​​WGS-84​​坐标系投影成正方形：

这是通过舍弃了南北​​85.051129纬度​​以上的地区实现的，因为它是正方形，所以一个大的正方形可以很方便的被分割为更小的正方形。

坐标系更详细的信息可参考​​GIS之坐标系统​​，​​EPSG:3857​​的详细信息可参考​​EPSG:3857​​。

经纬度定位行列号

上一节里我们简单介绍了一下坐标系，按照​​Web​​地图的标准，我们的地图引擎也选择支持​​EPSG:3857​​投影，但是我们通过高德工具获取到的是火星坐标系的经纬度坐标，所以第一步要把经纬度坐标转换为​​Web墨卡托​​投影坐标，这里为了简单，先直接把火星坐标当做​​WGS-84​​坐标，后面再来看这个问题。 转换方法网上一搜就有：

// 角度转弧度
const angleToRad = (angle) => 
    return angle * (Math.PI / 180)


// 弧度转角度
const radToAngle = (rad) => 
    return rad * (180 / Math.PI)


// 地球半径
const EARTH_RAD = 6378137

// 4326转3857
const lngLat2Mercator = (lng, lat) => 
    // 经度先转弧度，然后因为 弧度 = 弧长 / 半径 ，得到弧长为 弧长 = 弧度 * 半径 
    let x = angleToRad(lng) * EARTH_RAD; 
    // 纬度先转弧度
    let rad = angleToRad(lat)
    // 下面我就看不懂了，各位随意。。。
    let sin = Math.sin(rad)
    let y = EARTH_RAD / 2 * Math.log((1 + sin) / (1 - sin))
    return [x, y]


// 3857转4326
const mercatorTolnglat = (x, y) => 
    let lng = radToAngle(x) / EARTH_RAD
    let lat = radToAngle((2 * Math.atan(Math.exp(y / EARTH_RAD)) - (Math.PI / 2)))
    return [lng, lat]

​​3857​​坐标有了，它的单位是​​米​​，那么怎么转换成瓦片的行列号呢，这就涉及到​​分辨率​​的概念了，即地图上一像素代表实际多少米，分辨率如果能从地图厂商的文档里获取是最好的，如果找不到，也可以简单计算一下（如果使用计算出来的也不行，那就只能求助搜索引擎了），我们知道地球半径是​​6378137​​米，​​3857​​坐标系把地球当做正圆球体来处理，所以可以算出地球周长，投影是贴着地球赤道的：

所以投影成正方形的世界平面图后的边长代表的就是地球的周长，前面我们也知道了每一层级的瓦片数量的计算方式，而一张瓦片的大小一般是​​256*256​​像素，所以用地球周长除以展开后的世界平面图的边长就知道了地图上每像素代表实际多少米：

// 地球周长
const EARTH_PERIMETER = 2 * Math.PI * EARTH_RAD
// 瓦片像素
const TILE_SIZE = 256

// 获取某一层级下的分辨率
const getResolution = (n) => 
    const tileNums = Math.pow(2, n)
    const tileTotalPx = tileNums * TILE_SIZE
    return EARTH_PERIMETER / tileTotalPx

地球周长算出来是​​40075016.68557849​​，可以看到​​OpenLayers​​就是这么计算的：

​​3857​​坐标的单位是​​米​​，那么把坐标除以分辨率就可以得到对应的像素坐标，再除以​​256​​，就可以得到瓦片的行列号：

函数如下：

// 根据3857坐标及缩放层级计算瓦片行列号
const getTileRowAndCol = (x, y, z) => 
    let resolution = getResolution(z)
    let row = Math.floor(x / resolution / TILE_SIZE)
    let col = Math.floor(y / resolution / TILE_SIZE)
    return [row, col]

接下来我们把层级固定为​​17​​，那么分辨率​​resolution​​就是​​1.194328566955879​​，雷峰塔的经纬度转成​​3857​​的坐标为：​​[13374895.665697495, 3533278.205310311]​​，使用上面的函数计算出来行列号为：​​[43744, 11556]​​，我们把这几个数据代入瓦片的地址里进行访问：

https://webrd01.is.autonavi.com/appmaptile?x=43744&y=11556&z=17&lang=zh_cn&size=1&scale=1&style=8

一片空白，这是为啥呢，其实是因为原点不一样，​​4326​​和​​3857​​坐标系的原点在赤道和本初子午线相交点，非洲边上的海里，而瓦片的原点在左上角：

再来看下图会更容易理解：

​​3857​​坐标系的原点相当于在世界平面图的中间，向右为​​x​​轴正方向，向上为​​y​​轴正方向，而瓦片地图的原点在左上角，所以我们需要根据图上【绿色虚线】的距离计算出【橙色实线】的距离，这也很简单，水平坐标就是水平绿色虚线的长度加上世界平面图的一半，垂直坐标就是世界平面图的一半减去垂直绿色虚线的长度，世界平面图的一半也就是地球周长的一半，修改​​getTileRowAndCol​​函数：

const getTileRowAndCol = (x, y, z) => 
  x += EARTH_PERIMETER / 2     // ++
  y = EARTH_PERIMETER / 2 - y  // ++
  let resolution 以上是关于从零打造一个Web地图引擎的主要内容，如果未能解决你的问题，请参考以下文章
 
 从零打造一个Web地图引擎
 WebGL 渲染管线在 Web 3D 地图中的应用
 百度地图的瓦片规则参数
 Web（瓦片）地图的工作原理
 科普Web（瓦片）地图的工作原理
 Web版百度地图加载离线瓦片