⽡⽚切图算法以及并发切图实践（上篇）
互联⽹地图服务商的在线地图都通过⽡⽚的⽅式提供，称为⽡⽚地图服务。最常见的地图⽡⽚是图⽚格式的，现在有的地图服务商也提供了⽮量的⽡⽚数据，然后在⽤户端使⽤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 ，如像详细了解请参照原⽂。
具有唯⼀的⽡⽚等级（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）中。
坐标转换图解
从⾼德地图坐标转换图解中可以看出，⾼德地图的坐标转换具有以下特点：
所有坐标转换都在某⼀⽡⽚等级下进⾏，不同⽡⽚等级下的转换结果不同。
经纬度坐标可以直接转换为⽡⽚坐标和⽡⽚像素坐标。
待到重阳日还来就菊花的意思⽡⽚像素坐标需要结合其⽡⽚坐标才能得到该像素坐标的经纬度坐标。
坐标转换公式
参照此
经纬度坐标（lng, lat）转⽡⽚坐标（tileX， tileY）:
\[tileX=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\rfloor \]
\[tileY=\lfloor{(\frac{1}{2}-\frac{\ln(\tan(lat\times\pi/180)+\c(lat\times\pi/180))}{2\timesπ}})\times{2^{Level}}\rfloor \]
经纬度坐标（lng, lat）转像素坐标（pixelX, pixelY）
\[pixelX=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\times256\%256\rfloor \]
\[pixelY=\lfloor{(1-\frac{\ln(\tan(lat\times\pi/180)+\c(lat\times\pi/180))}{2\timesπ}})\times{2^{Level}}\times256\%256\rfloor \]
⽡⽚（tileX， tileY）的像素坐标（pixelX, pixelY）转经纬度坐标（lng, lat）
\[lng=\frac{tileX+\frac{pixelX}{256}}{2^{Level}}\times360-180 \]
\[lat=\arctan({\sinh({\pi-2\times\pi\times\frac{tileY+\frac{pixelY}{256}}{2^{Level}}})})\times\frac{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=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\rfloor \]
\[tileY=\lfloor{(\frac{1}{2}-\frac{\ln(\tan(lat\times\pi/180)+\c(lat\times\pi/180))}{2\timesπ}})\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，表⽰⼀个⽡⽚切图的长宽。
缩放图⽚
缩放图⽚我是⽤的是图⽚库。具体缩放代码如下。
tileImage = TileImage.Clone(c => { c.Resize(newSize); });
确定偏移
例如在上图中切图偏移量就是中间点距离左上⾓切图坐标的距离。也就是西北⾓的⽡⽚坐标的⼩数部分 * ⽡⽚切图的规格（长宽像素）
本手妙手俗手高考作文var xPixelOfft = -(int) ((minorPoint.X - Math.Truncate(minorPoint.X)) * TileExtensions.TileOpacity);
var yPixelOfft = -(int) ((minorPoint.Y - Math.Truncate(minorPoint.Y)) * TileExtensions.TileOpacity);
代码实现如上。
裁切
var (x, y) = tileMapCoordinate;
var xTilePixelOfft = x * TileExtensions.TileOpacity + xPixelOfft;
var yTilePixelOfft = y * TileExtensions.TileOpacity + yPixelOfft;
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 = xTilePixelOfft + pixelCol;
var pixelY = yTilePixelOfft + 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 表⽰此⽡⽚切图相对于西北⾓⽡⽚的相对坐标。在存储处理完的切⽚时需要额外加上西北⾓⽡⽚的绝对坐标，表⽰此⽡⽚的绝对坐标。贵州一本线
性能部分
性能部分将在下章描述。性能部分将着重描述并发，以及多线程粒度和内存管理。
