一、背景
最近工作上有這么一個需求：根據坐標點獲取Dem高程值。當前是通過Geoserver發布的WMTS柵格地圖服務，Geoserver示例中可以實現點擊地圖查詢點對應的像素值。如下圖所示，點擊地圖后左下角即可展示地圖點擊的柵格像素值。
二、資料查閱
既然示例代碼中已經有方法了，就直接F12查看請求的Url地址，在Url地址中發現了這四個比較關鍵的參數：TileCol=1735&TileRow=1036&I=36&J=252，示例請求地址如下圖：
接下來，去Geoserver官網用戶手冊查詢參數的意義，其中i和j是瓦片的像素坐標，TileCol和TileRow是對應的瓦片行列號。Geoserver根據經緯度查詢柵格值思路是先根據經緯度查找對應的瓦片行列號，再根據經緯度查找瓦片上的像素坐標，這樣就可以定位到那張瓦片上哪個像素點了，這些信息獲取到了自然就可以取到像素值了。以上就是大致的解決思路。
- Geoserver用戶手冊地址：https://www.osgeo.cn/geoserver-user-manual/services/wms/reference.html
- 論文GIS依據移步大佬地址：http://www.doc88.com/p-6116923730629.html
- 切片算法策略移步大佬博客：切圖算法簡述1——WGS 84坐標系下XYZ（WMTS）、TMS切片算法策略_金左手33的博客-CSDN博客
三、編碼驗證
根據經緯度和級別獲取瓦片行列號摘錄於大佬博客，了解更多可移步詳讀，關鍵代碼如下：
/// <summary> /// 獲取xyz(WMTS)服務瓦片左上角經緯度 /// </summary> /// <param name="column">列序號（x序列）</param> /// <param name="row">行序號（y序列）</param> /// <param name="zoom">地圖層級</param> /// <returns></returns> public static LatLngInfo XYZLatLng(int column, int row, int zoom) { double lon = XYZLongitude(column, zoom); double lat = XYZLatitude(row, zoom); return new LatLngInfo() { Lat = lat, Lng = lon }; } /// <summary> /// 通過列序列號、zoom層級獲取xyz(WMTS)服務瓦片左上角經度 /// </summary> /// <param name="column">列序列號（x序列）</param> /// <param name="zoom">地圖層級</param> /// <returns></returns> public static double XYZLongitude(int column, int zoom) { return (column * 180.0) / Math.Pow(2.0, zoom) - 180.0; } /// <summary> /// 通過行序列號、zoom層級獲取xyz(WMTS)服務瓦片左上角緯度 /// </summary> /// <param name="row">行序列號（y序列）</param> /// <param name="zoom">地圖層級</param> /// <returns></returns> public static double XYZLatitude(int row, int zoom) { double latitude = (row * 180.0) / Math.Pow(2.0, zoom) - 90.0; return -latitude; } /// <summary> /// 通過坐標經緯度、zoom層級獲取xyz(WMTS)服務瓦片行列號（xy） /// </summary> /// <param name="longitude">經度</param> /// <param name="latitude">緯度</param> /// <param name="zoom"></param> /// <returns></returns> public static TileIndexInfo XYZTileIndex(double longitude, double latitude, int zoom) { int column = XYZColumn(longitude, zoom); int row = XYZRow(latitude, zoom); return new TileIndexInfo() { Column = column, Row = row }; } /// <summary> /// 通過坐標經度、zoom層級獲取xyz(WMTS)服務瓦片列序列號（x序列） /// </summary> /// <param name="longitude">經度</param> /// <param name="zoom">地圖層級</param> /// <returns>列序列號（x序列）</returns> public static int XYZColumn(double longitude, int zoom) { return (int)Math.Floor((longitude + 180.0) / 180.0 * Math.Pow(2, zoom)); } /// <summary> /// 通過坐標緯度、zoom層級獲取xyz(WMTS)服務瓦片行序列號（y序列） /// </summary> /// <param name="latitude">緯度</param> /// <param name="zoom">地圖層級</param> /// <returns>行序列號（y序列）</returns> public static int XYZRow(double latitude, int zoom) { return (int)Math.Floor((-latitude + 90) / 180.0 * Math.Pow(2, zoom)); }
根據經緯度和級別獲取坐標點對應瓦片的像素坐標：
/// <summary> /// 計算坐標點對應瓦片上的像素坐標 /// </summary> /// <param name="coord1">當前坐標點</param> /// <param name="topLeftInfoCoord">瓦片左上角坐標</param> /// <param name="level">級別</param> /// <returns></returns> static (int i, int j) GetIJ(LatLngInfo coord1, LatLngInfo topLeftInfoCoord, int level) { var i = GetPixelX(coord1.Lng, level) - GetPixelX(topLeftInfoCoord.Lng, level); var j = GetPixelY(coord1.Lat, level) - GetPixelY(topLeftInfoCoord.Lat, level); return ((int)i, (int)j); } /// <summary> /// 根據經緯度獲取像素坐標 /// </summary> /// <param name="lon">經度</param> /// <param name="lat">維度</param> /// <param name="level">級別</param> /// <returns></returns> static (double px, double py) GetPixelXY(double lon, double lat, int level) { var px = Math.Round((lon + 180) / 360 * Math.Pow(2, level) * 256 % 256); var py = Math.Round((1 - Math.Log(Math.Tan(lat * Math.PI / 180) + 1 / Math.Cos(lat * Math.PI / 180)) / (2 * Math.PI)) * Math.Pow(2, level) * 256 % 256); return (px, py); }
四、總結
感謝以上大佬的知識分享，現簡單總結下通過經緯度獲取高程值的知識點：
- 先要通過經緯度獲取到瓦片的行列號，索引到具體的某張瓦片上；
- 再通過坐標點計算瓦片左上角的經緯度值，以備計算像素坐標用；
- 分別計算當前坐標點和瓦片左上角坐標點的像素坐標值；
- 兩個像素坐標相減即可得到瓦片上的像素坐標；
- 參數計算好拼接Url就可通過Http請求獲取到高程值。
參考資料：