TMS逆向到WMTS - 陈不易

有一個 TMS 的瓦片數據源，需要 “模擬” 一個 WMTS 服務出來，需要怎麼做？

這個情況，其實有現成的基礎設施或者說輪子來解決，比如各個地圖伺服器等，.net 生態也有 tile-map-service-net5這種開源工具，這個問題之所以是個問題在於兩個限制條件。

所用客戶端不支持加載 XYZ/TMS 格式的數據，只能加載 WMS 和 WMTS 格式的數據。
使用的數據是切好片的 TMS 結構的數據。
客戶端不方便依賴外部地圖伺服器。

模仿資源鏈接#

一些我們熟悉的互聯網地圖，用的都是 XYZ 或者 TMS 的方式，例如 OSM、Google Map 和 Mapbox 等等，從之前的栅格瓦片到如今矢量瓦片更為常見，想要用 TMS “模仿” WMTS 的請求格式，需要先了解他們直接有啥不一樣。

XYZ（slippy map tilename）#

256*256 像素的圖片
每個 Zoom 層級是一個文件夾，每個 Column 是個子文件夾，每個瓦片是一個用 Row 命名的圖片文件
格式類似/zoom/x/y.png
x 在 (180°W ~ 180°E)，y 在（85.0511°N ~85.0551°S），Y 軸從頂部向下。

可以從Openlayers TileDebug Example，看到一個簡單的 XYZ 瓦片的示例。

TMS#

TMS 的 Wiki wikipedia沒涉及什麼細節、osgeo-specification
只描述了協議的一些應用細節。反倒是 geoserver docs 關於 TMS 的部分寫的更務實一些。 TMS 是 WMTS 的前身，也是 OSGeo 制定的標準。

請求形如：
http://host-name/tms/1.0.0/layer-name/0/0/0.png

為了支持多種文件格式和空間參考系統，也可以指定多個參數：
http://host-name/tms/1.0.0/layer-name@griset-id@format-extension/z/x/y

TMS 標準的瓦片格網從左下角開始，Y 軸從底部向上。有的地圖伺服器，例如 geoserver，就支持一個額外的參數flipY=true 來翻轉 Y 坐標，這樣就可以兼容 Y 軸從頂部向下的服務類型，比如 WMTS 和 XYZ。

tms-grid

WMTS#

WMTS 相較上述兩個直觀的協議，內容更複雜，支持的場景也更多。2010 年由OGC第一次公布。起始在此之前，1997 年 Allan Doyle 的論文 “Www mapping framework” 之後，OGC 就開始謀劃網絡地圖相關標準的制定了。在 WMTS 之前，最早的，也是應用最廣泛的網絡地圖服務標準是 WMS。因為 WMS 每個請求是依據用戶地圖縮放級別和螢幕大小來組織地圖響應，這些響應大小各異，在多核 CPU 還沒那麼普及的當年，這種按需即時生成地圖的方式非常奢侈，同時想要提升響應速度非常困難。於是有開發者開始嘗試預先生成瓦片的方式，於是湧現出了許多方案，前面提到的 TMS 就是其中的一個，後面 WMTS 應運而生，開始被廣泛應用。 WMTS 支持鍵值對 (kvp) 和 Restful 的方式對請求參數編碼。

KVP 形如：
<baseUrl>/layer=<full layer name>&style={style}&tilematrixset={TileMatrixSet}}&Service=WMTS&Request=GetTile&Version=1.0.0&Format=<imageFormat>&TileMatrix={TileMatrix}&TileCol={TileCol}&TileRow={TileRow}
Restful 形如：
<baseUrl>/<full layer name>/{style}/{TileMatrixSet}/{TileMatrix}/{TileRow}/{TileCol}?format=<imageFormat>
由於是栅格瓦片，這裡只需要找到 XYZ 與瓦片矩陣和瓦片行列號的對應關係就好了。

TileMatrix
TileRow
TileCol

wmts-grid
這裡的瓦片行列號是從左上角開始的，Y 軸從頂部向下。

這樣，就找到了 TMS 與 WMTS 的各參數對應關係，接下來就是如何把 TMS 轉換成 WMTS 的請求了，如下：

TileRow = 2^zoom - 1 - y = (1 << zoom) - 1 - y
TileCol = x
TileMatrix = zoom

在不考慮其他空間參考的情況下，縮放層級對應瓦片矩陣，x 對應瓦片列號，y 取反（因為起始方向相反）。

模擬一個 WMTS Capabilities 描述文件#

WMTS 規範的要求，幾乎可以說是細到頭髮絲，所以各個客戶端，不管是 Web 端的 Openlayers ，還是桌面端的 QGIS 或 Skyline 等，都支持直接解析 Capabilities 描述文件，然後根據描述文件的內容來選擇圖層、樣式和空間參考，所以我們這裡還要模擬一個 WMTS Capabilities 描述文件出來。

Capabilities 描述文件的構成#

一個 WMTS Capabilities 描述文件的例子可以在opengis schema,天地圖山東找到。

Capabilities 描述文件的內容非常多，這裡只列出一些重要的部分（忽略標題，聯繫方式等）：

OperationsMetadata:
  - GetCapabilities >> 獲取 Capabilities 描述文件的方式
  - GetTile >> 獲取瓦片的方式

Contents:
  - Layer
    - boundingBox >> 圖層的經緯度範圍
    - Style
    - TileMatrixSetLink >> 圖層支持的空間參考
    - TileMatrixSet >> 空間參考
      - TileMatrixSetLimits >> 空間參考的縮放層級範圍
        - TileMatrixLimits >> 每個縮放層級的瓦片行列號範圍
  - Style
  - TileMatrixSet
    - TileMatrix

關鍵的部分就是 boundingBox、TileMatrixSetLimits、TileMatrixLimits ，只需要根據圖層的空間參考和縮放層級來計算出來就好了。

boundingBox 的計算比較簡單，就是圖層的經緯度範圍，這裡就不展開了。

TileMatrixSetLimits 的計算比較簡單，就是圖層的空間參考的縮放層級範圍。

TileMatrixLimits 的計算比較複雜，可以只在圖層範圍比較小的時候再弄，全球地圖就沒必要了，需要根據圖層的空間參考和縮放層級來計算出來，下面是一段偽代碼（4326 到 3857）。

FUNCTION GetTileRange(minLon, maxLon, minLat, maxLat, zoom, tile_size = 256)

minLonRad = minLon * PI / 180
maxLonRad = maxLon * PI / 180
minLatRad = minLat * PI / 180
maxLatRad = maxLat * PI / 180

tile_min_x = Floor((minLonRad + PI) / (2 * PI) * Pow(2, zoom))
tile_max_x = Floor((maxLonRad + PI) / (2 * PI) * Pow(2, zoom))
tile_min_y = Floor((PI - Log(Tan(minLatRad) + 1 / Cos(minLatRad))) / (2 * PI) * Pow(2, zoom))
tile_max_y = Floor((PI - Log(Tan(maxLatRad) + 1 / Cos(maxLatRad))) / (2 * PI) * Pow(2, zoom))

// adjust tile range based on tile size
tile_min_x = Floor((double)tile_min_x * tile_size / 256)
tile_max_x = Ceiling((double)tile_max_x * tile_size / 256)
tile_min_y = Floor((double)tile_min_y * tile_size / 256)
tile_max_y = Ceiling((double)tile_max_y * tile_size / 256)

RETURN  (tile_min_x, tile_max_x, tile_min_y, tile_max_y)

生成 WMTS Capabilities 描述文件#

生成一個最小化的 WMTS Capabilities 描述文件，把上面的關鍵部分填充上，之後構造一個指向標準描述文件地址的的 Restful 風格的 URL。

後話#

以上是一個簡單的 TMS 轉 WMTS 的思路，實際上還有很多細節需要考慮，比如空間參考的轉換，縮放層級的轉換，瓦片行列號的轉換，瓦片的格式轉換等等。
期間也踩了一些坑，感覺這部分更有意思。

第一部分，很快就參考 tile-map-service-net5 的思路，完成了 y >> tileRow的轉換。代碼在WebMercator.cs ，其實在 StackOverflow 上也有人問過這個問題，是有答案的，但我還是選擇從軟件裡找答案，因為這樣自己心裡更踏實。

第二部分就很頭大，首先模擬出了資源鏈接，構建了一個簡單的 XML，但是在目標客戶端上不能直接加載，很直接的想到了通過標準服務測試一下，然後哪來一個 Capabilities 描述文件來修改。己想首先在比較熟悉的 Openlayers 上測試，然後再去修改 Capabilities 描述文件。Openlayers 的加載方式還是很靈活的，在沒有 Capabilities 描述文件的情況下，可以直接通過配置參數訪問。

// fetch the WMTS Capabilities parse to the capabilities 
const options = optionsFromCapabilities(capabilities, {
            layer: 'nurc:Pk50095',
            matrixSet: 'EPSG:900913',
            format: 'image/png',
            style: 'default',
});
const wmts_layer =new TileLayer({
    opacity: 1,
    source: new  WMTS(options),
})

很遺憾，瓦片沒有加載上，甚至networks裡沒有發送請求。於是又去另一個 WMTS 相關的例子哪裡，自定義了一個 TileGrid，然後把瓦片的行列號轉換成了 3857 的行列號，這時候可以加載了。

const projection = getProjection('EPSG:3857');
const projectionExtent = projection.getExtent();
const size = getWidth(projectionExtent) / 256;
const resolutions = new Array(31);
const matrixIds = new Array(31);
for (let z = 0; z < 31; ++z) {
    // generate resolutions and matrixIds arrays for this WMTS
    resolutions[z] = size / Math.pow(2, z);
    matrixIds[z] = `EPSG:900913:${z}`;
}
var wmtsTileGrid = new WMTSTileGrid({
    origin: getTopLeft(projectionExtent), resolutions: resolutions, matrixIds: matrixIds,
})

在確認了是 TileGrid 的問題之後，首先將自己生成的 TileGrid 與 Openlayers 從 Capabilities 解析出來的 TileGrid 進行對比。發現自己生成的 TileGrid 有一些字段是空的，於是挨個測試，最後發現設置fullTileRanges_和extent_兩個內部參數為空時，影像可以加載。

去翻 OL 源碼，發現fullTileRanges_和extent_在 getFullTileRange中被用到。

也就是說，當fullTileRanges_和extent_為空時，getFullTileRange會返回一個空的範圍。

而getFullTileRange在withinExtentAndZ中用到了，這裡是用來判斷當前可視區域是否有該圖層的瓦片。
也就是說，當fullTileRanges_和extent_為空時，獲取不到 TileRange，withinExtentAndZ會一直返回true，這樣就會一直加載瓦片了，也就是加載成功的原因。

相反，從 Capabilities 解析出來的fullTileRanges_和extent_指向了錯誤的TileRange，導致withinExtentAndZ一直返回false，這樣就不會加載瓦片了，也就是加載失敗的原因。

終於找到了原因，但這裡又被騙了。在wmts.js，構造函數上有一行註釋:

class WMTS extends TileImage {
  /**
   * @param {Options} options WMTS options.
   */
  constructor(options) {
      // TODO: add support for TileMatrixLimits
  }
}

這使我開始的時候誤以為，fullTileRanges_和extent_是根據經緯度範圍（boundingBox）計算出來的，而不是根據TileMatrixLimits算的，於是乎又檢查了一遍 boundingBox，確認無誤後，才開始著手修改TileMatrixLimits。

開始的時候，以為 TileMatrixLimits 是每個層級的瓦片範圍，而不是圖層的範圍，所以沒注意到這個參數，這才走了彎路。

寫在 2023 年，WMTS 已經不是一個新的協議了，OGC Tile API 已經成為正式標準了，自己對 WMTS 了解還是半瓶水，真是汗顏😅。