瓦片切圖算法以及並發切圖實踐(上篇)

互聯網地圖服務商的在線地圖都通過瓦片的方式提供，稱為瓦片地圖服務。最常見的地圖瓦片是圖片格式的，現在有的地圖服務商也提供了矢量的瓦片數據，然后在用戶端使用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=\lfloor\frac{lng + 180}{360}\times{2^{Level}}\rfloor \]

\[tileY=\lfloor{(\frac{1}{2}-\frac{\ln(\tan(lat\times\pi/180)+\sec(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)+\sec(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)+\sec(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，表示一個瓦片切圖的長寬。

縮放圖片

縮放圖片我是用的是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 表示此瓦片切圖相對於西北角瓦片的相對坐標。在存儲處理完的切片時需要額外加上西北角瓦片的絕對坐標，表示此瓦片的絕對坐標。

性能部分

性能部分將在下章描述。性能部分將着重描述並發，以及多線程粒度和內存管理。