一、背景
最近工作上有这么一个需求：根据坐标点获取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请求获取到高程值。
参考资料：