瓦片地图服务与地图瓦片原理
搜
刘艳鹏约 14 分钟
刘艳鹏约 14 分钟一、瓦片地图服务
二、地图瓦片原理
2.1 瓦片
2.2 金字塔模型
根据WMTS官方标准，规定的瓦片地图设定Web-Mercator投影坐标系原点为东经180°、南纬85.05°，x轴向右，y轴向上，zoom level最小为0，最大为24。瓦片尺寸通常为256*256像素（也可以是其它尺寸）。
但是在这个标准出台以前，Google的瓦片地图规范已经施行多年，深受行业认可，国内外图商，如OSM、高德、腾讯等都使用Google的瓦片地图规范，其定义略有不同，坐标原点为东经180°、北纬85.05°，x轴向右，y轴向下，zoom level最小为0、最大为21。
三、瓦片地图的生产与分发
3.1 地图投影与配图
3.1.1 数据投影
（1）地图坐标系
地理坐标系，也叫经纬度坐标系。像素坐标系，也可以称为屏幕坐标系。像素坐标系和地理坐标系存在对应关系，屏幕上的每一个像素都对应一个经纬度位置。在瓦片地图服务中，不同的zoom level下，像素坐标系与地理坐标系之间的对应关系是不同的。 投影坐标系，其定义在上面已经给出，既然经纬度能够在二维平面进行表示，那么，就可以在二维平面上对像素坐标系与投影坐标系之间进行映射，这就引出了投影坐标系的一个重要作用：将地理坐标系和屏幕坐标系关联起来。所以，理解投影坐标系是GIS研究中的一个重要环节。
（2）投影方法与分类
（3）墨卡托投影
图9：手机与地图的位置关系
3.1.2 地图配图
3.2 地图切片及瓦片定位
下面我们以Google瓦片规范简单描述一下Web墨卡托投影的算法原理。当缩放等级zoom level等于0时，Web墨卡托投影就是球面投影到一个正方形的平面上，这个平面就是将世界坐标调整为左上角为(0,0)，右下角为(256,256)的正方形，假设基本单位为像素，那么就等于把地图投影在一个256pixel*256pixel的图幅上：
x = \lfloor \frac{256}{2\pi} * 2^{zoom\_level}(lng + \pi) \rfloor * pixels
\\
y = \lfloor \frac{256}{2\pi} * 2^{zoom\_level}(\pi - \ln[\tan(\frac{\pi}{4} + \frac{lat}{2})]) * pixels
3.2.1 地理坐标=>投影坐标（经纬度=>米）
X_{proj} = lng * \frac{2 * \pi * R}{2.0}
\\
Y_{proj} = R * log(tan(\frac{(90+lat)*\pi}{360}))
# 地理坐标投影
def LatLongToMeterXY(lng, lat) {
# 地球的周长的一半 20037508.342789244 单位米
circumferenceHalf = math.pi * 2 * 6378137 / 2.0
meterX = lng * circumferenceHalf / 180
temp = math.log(math.tan((90 + lat) * (math.pi / 360.0))) / (math.pi / 180.0)
meterY = circumferenceHalf * temp / 180
return meterX,meterY
# [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
3.2.2 投影坐标=>像素坐标（米=>像素）
X_{pixel} = \frac{(X_{proj}+\pi*R)*tile_size*z^{zoom}}{2*\pi*R}
\\
Y_{pixel} = \frac{(Y_{proj}+\pi*R)*tile_size*z^{zoom}}{2*\pi*R}
# 投影坐标转像素坐标
def metersToPixelXY(meterX, meterY, zoom, tileSize=256) {
# 地球的周长的一半 20037508.342789244 单位米
circumferenceHalf = math.pi * 2 * 6378137 / 2.0
# 米/每像素
resolution = math.pi * 2 * 6378137 / (tileSize * math.pow(2, zoom))
pixelX = (meterX + circumferenceHalf) / resolution
pixelY = (meterY + circumferenceHalf) / resolution
return pixelX,pixelY
3.2.3 地理坐标=>像素坐标（经纬度=>像素）
X_{pixel} = \frac{lng+180}{360}*2^z*256 \mod 256
\\
Y_{pixel} = (1-\frac{\ln(\tan(\frac{lat*\pi}{180}) + \sec(\frac{lat * \pi}{180}))}{2* \pi} *2^z*256\mod256
# 经纬度转像素
def lnglatToPixel(lng,lat,zoom):
pixelX=round((lng+180)/360*math.pow(2,zoom)*256%256)
pixelY=round((1-math.log(math.tan(math.radians(lat))+1/math.cos(math.radians(lat)))/(2*math.pi))*math.pow(2,zoom)*256%256)
return pixelX,pixelY
3.2.4 像素坐标 =>瓦片坐标（像素=>瓦片）
X_{tile} = \lfloor \frac{X_{pixel}}{tile\_size} \rfloor
\\
Y_{tile} = \lfloor \frac{Y_{pixel}}{tile\_size} \rfloor
def pixelXYToTileXY(pixelX, pixelY, tileSize=256) {
tileX = Math.floor(pixelX / tileSize)
tileY = Math.floor(pixelY / tileSize)
return tileX,tileY
3.2.5 地理坐标<=>瓦片坐标（经纬度<=>瓦片）
（1）地理坐标转瓦片坐标
X_{pixel} = \frac{lng+180}{360}*2^z*256 \mod 256
\\
Y_{pixel} = (1-\frac{\ln(\tan(\frac{lat*\pi}{180}) + \sec(\frac{lat * \pi}{180}))}{2* \pi} *2^z*256\mod256
# 经纬度转瓦片
def lnglatToTile(lng,lat,zoom):
tileX=int((lng+180)/360*math.pow(2,zoom))
tileY=int((1-math.asinh(math.tan(math.radians(lat)))/math.pi)*math.pow(2,zoom-1))
return tileX,tileY
（2）瓦片坐标转地理坐标
lng = \frac{X_{tile}}{2^z}*360-180
\\
lat = arctan(sinh(\pi*(1-\frac{2*Y_{tile}}{2^z})))*\frac{180}{\pi}
# 瓦片转经纬度
def tileToLnglat(tileX,tileY,zoom):
lng=tileX/math.pow(2,zoom)*360-180
lat=math.degrees(math.atan(math.sinh(math.pi*(1-2*tileY/math.pow(2,zoom)))))
return lng,lat
（3）瓦片坐标的像素坐标转地理坐标
lng = \frac{X_{tile}+\frac{X_{pixel}}{256}}{2^z}*360-180
\\
lat = arctan(sinh(\pi*(1-\frac{2*(Y_{tile}+{\frac{Y_{pixel}}{256}})}{2^z})))*\frac{180}{\pi}
# 瓦片坐标的像素坐标转经纬度
def pixelToLnglat(tileX,tileY,pixelX,pixelY,level):
lng=(tileX+pixelX/256)/math.pow(2,level)*360-180
lat=math.degrees(math.atan(math.sinh(math.pi-2*math.pi*(tileY+pixelY/256)/math.pow(2,level))))
return lng,lat
3.3 瓦片分发
A 子域名：https://a.tile.openstreetmap.org/12/3371/1551.png B 子域名：https://b.tile.openstreetmap.org/12/3371/1551.png C 子域名：https://c.tile.openstreetmap.org/12/3371/1551.png
通过上面的URL以及瓦片路径结构，可以看出上图中的瓦片来自于缩放等级为12，第3771列、第1551行。可以概括的看这个URL结构为：
https://{sub_domain}.tile.openstreetmap.org/{z}/{x}/{y}.png
3.4 瓦片地图服务解析
3.4.1 瓦片地图服务解析流程
3.4.2 leaflet中实现瓦片的渲染
<html>
<head>
<title>瓦片原理-OSM</title>
<meta charset="utf-8" />
<link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css" integrity="sha256-kLaT2GOSpHechhsozzB+flnD+zUyjE2LlfWPgU04xyI=" crossorigin="" />
<script src="https://unpkg.com/[email protected]/dist/leaflet.js" integrity="sha256-WBkoXOwTeyKclOHuWtc+i2uENFpDZ9YPdf5Hf+D7ewM=" crossorigin=""></script>
<style>
body {
padding: 0;
margin: 0;
}
html, body, #map {
height: 100%;
width: 100%;
}
</style>
</head>
<body>
<div id="map"></div>
</body>
<script>
var map = new L.Map('map', { center: [40, 116], zoom: 17});
var tiles = new L.GridLayer();
tiles.createTile = function(coords) {
console.log('coords:');
console.log(coords);
// 创建DOM装载瓦片
var tile = L.DomUtil.create('canvas', 'leaflet-tile');
var ctx = tile.getContext('2d');
// 获取瓦片尺寸
var size = this.getTileSize();
console.log('size:');
console.log(size);
tile.width = size.x;
tile.height = size.y;
// 计算瓦片左上角的像素坐标
var nwPoint = coords.scaleBy(size);
// 计算瓦片左上角的地理坐标
var nw = map.unproject(nwPoint, coords.z);
// 计算瓦片右下角的像素坐标
seCoords = new L.Point(coords.x+1,coords.y+1);
seCoords.z = coords.z;
console.log('seCoords:');
console.log(seCoords);
// 计算瓦片右下角的地理坐标
var sePoint = seCoords.scaleBy(size);
var se = map.unproject(sePoint, coords.z);
// 加载瓦片图像
var img = new Image();
var subdomains = ['a', 'b', 'c'];
var index = Math.abs(coords.x + coords.y) % subdomains.length;
img.src = 'http://'+subdomains[index]+'.tile.openstreetmap.org/'+coords.z+'/'+coords.x+'/'+coords.y+'.png';
console.log(img);
ctx.drawImage(img,0,0,size.x,size.y);
console.log(ctx);
// 标识对应的瓦片信息
// 背景颜色
ctx.fillStyle = 'white';
// 背景尺寸
ctx.fillRect(0, 0, size.x, 90);
// 字体颜色
ctx.fillStyle = 'black';
ctx.fillText('x: ' + coords.x + ', y: ' + coords.y + ', zoom level: ' + coords.z, 20, 20);
ctx.fillText('minlng: ' + nw.lng + ', minlat: ' + se.lat, 20, 40);
ctx.fillText('maxlng: ' + se.lng + ', matlat: ' + nw.lat, 20, 60);
var dlng = se.lng - nw.lng;
var dlat = nw.lat - se.lat;
ctx.fillText('dlng: ' + dlng + ', dlat: ' + dlat, 20, 80);
ctx.strokeStyle = 'red';
ctx.beginPath();
ctx.moveTo(0, 0);
ctx.lineTo(size.x-1, 0);
ctx.lineTo(size.x-1, size.y-1);
ctx.lineTo(0, size.y-1);
ctx.closePath();
ctx.stroke();
return tile;
};
// 加载OSM瓦片代码
// L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
// attribution: 'Map data © <a href="http://www.osm.org">OpenStreetMap</a>'
// }).addTo(map);
tiles.addTo(map);
</script>
</html>
四、总结
参考文献
[2]https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
[3]https://www.yuque.com/antv/l7/up4re5
[4]http://webgis.cn/standard-wmts.html
[5]https://github.com/antvis/L7/tree/master/packages/utils/src/tileset-manager