瓦片切图算法以及并发切图实践(上篇)

  • A+
所属分类：.NET技术
摘要

互联网地图服务商的在线地图都通过瓦片的方式提供，称为瓦片地图服务。最常见的地图瓦片是图片格式的，现在有的地图服务商也提供了矢量的瓦片数据，然后在用户端使用Canvas渲染成图片，如node-canvas实现百度地图个性化底图绘制。
在进行地图开发时，为获取特定经纬度所在区域的瓦片和获取瓦片上像素点对应的经纬度，经常需要进行经纬度坐标与瓦片坐标、像素坐标的相互转换。
便宜好用的国外VPS推荐

互联网地图服务商的在线地图都通过瓦片的方式提供，称为瓦片地图服务。最常见的地图瓦片是图片格式的，现在有的地图服务商也提供了矢量的瓦片数据，然后在用户端使用Canvas渲染成图片，如node-canvas实现百度地图个性化底图绘制。
在进行地图开发时，为获取特定经纬度所在区域的瓦片和获取瓦片上像素点对应的经纬度，经常需要进行经纬度坐标与瓦片坐标、像素坐标的相互转换。

基础知识

瓦片地图?️

在Gis地图开发中，经常会遇到超过屏幕大小的地图。 例如无人机给了一张巨大的航拍图片，如何将这张图片在不同的地图API中访问（高德地图、OpenLayers等等)。

瓦片地图就是为了解决此类问题产生的，一张大的世界地图或者背景图可以由几种地形来表示，每种地形对应一张小的的图片，我们称这些小的地形图片为瓦片。把这些瓦片拼接在一起，一个完整的地图就组合出来了，这就是瓦片地图的原理。在高德地图中每次放大缩小滚轮中，其实会加载的每一层级的TileMaps（瓦片切图），而我们需要做的就是需要将一张巨大的航拍地图，处理成一层一层的瓦片切图。

瓦片切图算法以及并发切图实践(上篇)

瓦片地图金字塔模型

瓦片地图金字塔模型是一种多分辨率层次模型，从瓦片金字塔的底层到顶层，分辨率越来越低，但表示的地理范围不变。 首先确定地图服务平台所要提供的缩放级别的数量 N，把缩放级别最高、地图比例尺最大的地图图片作为金字塔的底层，即第 0 层，并对其进行分块，从地图图片的左上角开始，从左至右、从上到下进行切割，分割成相同大小的正方形地图瓦片，形成第 0 层瓦片矩阵；在第 0 层地图图片的基础上，按每像素分割为 2×2 个像素的方法生成第 1 层地图图片，并对其进行分块，分割成与下一层相同大小的正方形地图瓦片，形成第1层瓦片矩阵；采用同样的方法生成第 2 层瓦片矩阵；...... 如此下去，直到第 N-1 层，构成整个瓦片金字塔。

瓦片切图算法以及并发切图实践(上篇)

瓦片地图一般采用ZXY规范的地图瓦片。（瓦片层级、瓦片 x坐标、瓦片 y坐标）

墨卡托投影

使用经纬度表示位置的大地坐标系虽然可以描述地球上点的位置，但是对于地图地理数据在二维平面内展示的场景，需要通过投影的方式将三维空间中的点映射到二维空间中。地图投影需要建立地球表面点与投影平面点的一一对应关系，在互联网地图中常使用墨卡托投影。墨卡托投影是荷兰地理学家墨卡托于1569年提出的一种地球投影方法，该方法是圆柱投影的一种。投影的更多内容，可以查看地图投影的N种姿势。

瓦片切图算法以及并发切图实践(上篇)

注意⚠️

  • 墨卡托投影并不是一种坐标系，而是为了在二维平面上展示三维地球而进行的一种空间映射。所以在GIS地图和互联网地图中，虽然用户看到的地图经过了墨卡托投影，但依然使用经纬度坐标来表示地球上点的位置。
    在地图绘制和地图可视化时，就需要将地图数据使用投影的方式来呈现。
  • 各大地图服务商都采用了Web Mercator进行投影，瓦片坐标系的不同主要是投影截取的地球范围不同、瓦片坐标起点不同。

瓦片切图规范

本文使用的是国际标准的经纬度坐标WGS84，Open Street Map。
此部分将大量翻译与参照此Wiki Slippy map tilenames，如像详细了解请参照原文。

  • 具有唯一的瓦片等级（Level）和瓦片坐标编号（tileX, tileY）。
  • 瓦片切图应当是 256*256 像素的一批PNG文件。
  • 每个缩放级别都是目录，每列是一个子目录，并且该列中的每个图块都是文件。
  • 瓦片切图的文件名（Url）的格式是 /zoom/x/y.png

瓦片地图等级范围

瓦片地图等级范围反映了地图可缩放的程度。
虽然最小的瓦片等级是0，但是部分地图并不提供0级或其他较小瓦片等级的地图，因为此时的世界地图将会很小，不能铺满用户设备窗口。

高德图片瓦片的层级是[1~19]

高德地图瓦片坐标

高德地图瓦片坐标与Google Map、Open Street Map相同。高德地图的墨卡托投影截取了纬度（约85.05ºS, 约85.05ºN）之间部分的地球，使得投影后的平面地图水平方向和垂直方向长度相等。将墨卡托投影地图的左上角作为瓦片坐标系起点，往左方向为X轴，X轴与北纬85.05º重合且方向向左；往下方向为Y轴，Y轴与东经180º（亦为西经180º）重合且方向向下。瓦片坐标最小等级为0级，此时平面地图是一个像素为256*256的瓦片。在某一瓦片层级Level下，瓦片坐标的X轴和Y轴各有(2^{Level})个瓦片编号，瓦片地图上的瓦片总数为(2^{Level}times2^{Level})。

瓦片切图算法以及并发切图实践(上篇)

如上图所示，此时X方向和Y方向各有4个瓦片编号，总瓦片数为16，层级2。中国大概位于高德瓦片坐标的（3，1）中。

坐标转换图解

瓦片切图算法以及并发切图实践(上篇)

从高德地图坐标转换图解中可以看出，高德地图的坐标转换具有以下特点：

  • 所有坐标转换都在某一瓦片等级下进行，不同瓦片等级下的转换结果不同。
  • 经纬度坐标可以直接转换为瓦片坐标和瓦片像素坐标。
  • 瓦片像素坐标需要结合其瓦片坐标才能得到该像素坐标的经纬度坐标。

坐标转换公式

参照此Slippy map tilenames

  • 经纬度坐标（lng, lat）转瓦片坐标（tileX， tileY）:
[tileX=lfloorfrac{lng + 180}{360}times{2^{Level}}rfloor ]

[tileY=lfloor{(frac{1}{2}-frac{ln(tan(lattimespi/180)+sec(lattimespi/180))}{2timesπ}})times{2^{Level}}rfloor ]

  • 经纬度坐标（lng, lat）转像素坐标（pixelX, pixelY）
[pixelX=lfloorfrac{lng + 180}{360}times{2^{Level}}times256%256rfloor ]

[pixelY=lfloor{(1-frac{ln(tan(lattimespi/180)+sec(lattimespi/180))}{2timesπ}})times{2^{Level}}times256%256rfloor ]

  • 瓦片（tileX， tileY）的像素坐标（pixelX, pixelY）转经纬度坐标（lng, lat）
[lng=frac{tileX+frac{pixelX}{256}}{2^{Level}}times360-180 ]

[lat=arctan({sinh({pi-2timespitimesfrac{tileY+frac{pixelY}{256}}{2^{Level}}})})timesfrac{180}{pi} ]

切图实践

有了上述的这些知识便可以编写我们切图的算法代码了。当业务中给我们一个高清的航拍图以及其西北角和东南角坐标时，我们需要如何通代码处理成瓦片切图。

瓦片切图算法以及并发切图实践(上篇)

端点到瓦片切图坐标的转换

var p = new Point {     X = (lngLat.Lng + 180.0) / 360.0 * (1 << zoom),     Y = (1.0 - Math.Log(Math.Tan(lngLat.Lat * Math.PI / 180.0) +                         1.0 / Math.Cos(lngLat.Lat * Math.PI / 180.0)) / Math.PI) / 2.0 * (1 << zoom) };  return p;

[tileX=lfloorfrac{lng + 180}{360}times{2^{Level}}rfloor ]

[tileY=lfloor{(frac{1}{2}-frac{ln(tan(lattimespi/180)+sec(lattimespi/180))}{2timesπ}})times{2^{Level}}rfloor ]

这部分代码是上述公式由转化过来的。但笔者代码中并未向下取整，是因为笔者想取得此经纬度转换到瓦片坐标的精确点位。比如当使用此算法处理西北角坐标时，会得出一个瓦片坐标，其x、y可能都是小数。

瓦片切图算法以及并发切图实践(上篇)

用瓦片切图来演示是如上图的效果，上图中间的点表示的则是航拍图的西北角。而此瓦片切图中超出西北角的部分为白色。对应到坐标系，因为将墨卡托投影地图的左上角作为瓦片坐标系起点。所以瓦片切图左上角的坐标，则代表此瓦片切图的坐标。

确定当前层级下所切图的分辨率

已经分别取得了西北和东南角的瓦片切图坐标。下一步便是要取得当前层级的所切图片。

Width = (int) Math.Floor((EastSouth.X - WestNorth.X) * TileExtensions.TileOpacity), Height = (int) Math.Floor((EastSouth.Y - WestNorth.Y) * TileExtensions.TileOpacity)

TileOpacity 默认为256，表示一个瓦片切图的长宽。

缩放图片

缩放图片我是用的是SixLabors.ImageSharp图片库。具体缩放代码如下。

tileImage = TileImage.Clone(c => { c.Resize(newSize); });

确定偏移

瓦片切图算法以及并发切图实践(上篇)

例如在上图中切图偏移量就是中间点距离左上角切图坐标的距离。也就是西北角的瓦片坐标的小数部分 * 瓦片切图的规格（长宽像素）

var xPixelOffset = -(int) ((minorPoint.X - Math.Truncate(minorPoint.X)) * TileExtensions.TileOpacity); var yPixelOffset = -(int) ((minorPoint.Y - Math.Truncate(minorPoint.Y)) * TileExtensions.TileOpacity);

代码实现如上。

裁切

var (x, y) = tileMapCoordinate;     var xTilePixelOffset = x * TileExtensions.TileOpacity + xPixelOffset; var yTilePixelOffset = y * TileExtensions.TileOpacity + yPixelOffset;  using var destImage = new Image<T>(TileExtensions.TileOpacity, TileExtensions.TileOpacity);  for (int pixelCol = 0; pixelCol < TileExtensions.TileOpacity; pixelCol++) {     for (int pixelRow = 0; pixelRow < TileExtensions.TileOpacity; pixelRow++)     {         var pixelX = xTilePixelOffset + pixelCol;         var pixelY = yTilePixelOffset + pixelRow;         tileImage.CopyTo(destImage, (pixelX, pixelY), (pixelCol, pixelRow));     } }  if (await TileImageRepository.TrySaveAsync(         new TileImage<T>(             tileId, destImage, tileZoomLevel,             (int) tilePositionInVertex.WestNorth.X + x,             (int) tilePositionInVertex.WestNorth.Y + y))) {     Interlocked.Increment(ref _finishCount); }

上述代码为单个瓦片切图的代码。x，y 表示此瓦片切图相对于西北角瓦片的相对坐标。在存储处理完的切片时需要额外加上西北角瓦片的绝对坐标，表示此瓦片的绝对坐标。

性能部分

性能部分将在下章描述。性能部分将着重描述并发，以及多线程粒度和内存管理。

便宜好用的国外VPS推荐