切图算法简述3——Web墨卡托（WebMercator）坐标系下XYZ（WMTS）、TMS切片算法策略

• 简述

 基于Web墨卡托平面坐标系就WMTS（XYZ）服务、TMS服务切图算法进行简述，下属内容来自于https://github.com/AliFlux/VectorTileRenderer/blob/master/VectorTileRenderer/GlobalMercator.cs中的算法整理，本人仅仅是做了总结，并且简单以C# 代码进行实现。各语言版本的开源地址连接在这儿：https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/，具体截图如下。

• 具体代码实现（C#）

 

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AlgorithmUtil.SRS
{
    /// <summary>
    /// 平面坐标实体类
    /// </summary>
    public class CoordinateInfo
    {
        public double X { get; set; }
        public double Y { get; set; }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AlgorithmUtil.SRS
{
    /// <summary>
    /// 经纬度坐标实体类
    /// </summary>
    public class LatLngInfo
    {
        public double Lat { get; set; }
        public double Lng { get; set; }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AlgorithmUtil.SRS
{
    /// <summary>
    /// 瓦片行列号实体类
    /// </summary>
    public class TileIndexInfo
    {
        public int Column { get; set; }
        public int Row { get; set; }

        public override string ToString()
        {
            return string.Format("{0},{1}", Column, Row);
        }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AlgorithmUtil.SRS
{
    public class GeoExtent
    {
        public double North { get; set; }
        public double South { get; set; }
        public double East { get; set; }
        public double West { get; set; }
    }
}

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AlgorithmUtil.SRS
{
    public class WebMercator
    {
        //瓦片大小
        private int tileSize;
        //初始分辨率
        private double initialResolution;
        private double originShift;

        public WebMercator()
        {
            this.tileSize = 256;
            this.initialResolution = 2 * Math.PI * 6378137 / tileSize;
            this.originShift = 2 * Math.PI * 6378137 / 2.0;
        }

        /// <summary>
        /// 根据地图层级获取地图分辨率
        /// </summary>
        /// <param name="zoom"></param>
        /// <returns></returns>
        private double Resolution(int zoom)
        {
            return this.initialResolution / (1 << zoom);
        }

        /// <summary>
        /// WGS84地理坐标转web墨卡托平面坐标
        /// </summary>
        /// <param name="lat">纬度</param>
        /// <param name="lon">经度</param>
        /// <returns>web墨卡托坐标</returns>
        public CoordinateInfo LatLonToMeters(double lat, double lon)
        {
            CoordinateInfo retval = new CoordinateInfo();
            try
            {
                retval.X = lon * this.originShift / 180.0;
                retval.Y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360.0)) / (Math.PI / 180.0);

                retval.Y *= this.originShift / 180.0;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// web墨卡托平面坐标转WGS84地理坐标
        /// </summary>
        /// <param name="mx">平面坐标x（单位：meter）</param>
        /// <param name="my">平面坐标y（单位：meter）</param>
        /// <returns>WGS84地理坐标</returns>
        public LatLngInfo MetersToLatLon(double mx, double my)
        {
            LatLngInfo latLng = new LatLngInfo();
            try
            {
                latLng.Lng = (mx / this.originShift) * 180.0;
                latLng.Lat = (my / this.originShift) * 180.0;

                latLng.Lat = 180 / Math.PI * (2 * Math.Atan(Math.Exp(latLng.Lat * Math.PI / 180.0)) - Math.PI / 2.0);
                return latLng;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 地图像素坐标转web墨卡托坐标
        /// </summary>
        /// <param name="px">像素坐标x</param>
        /// <param name="py">像素坐标y</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns>web墨卡托坐标</returns>
        public CoordinateInfo PixelsToMeters(double px, double py, int zoom)
        {
            CoordinateInfo retval = new CoordinateInfo();
            try
            {
                var res = Resolution(zoom);
                retval.X = px * res - this.originShift;
                retval.Y = py * res - this.originShift;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// web墨卡托坐标转地图像素坐标
        /// </summary>
        /// <param name="mx">平面坐标x（单位：meter）</param>
        /// <param name="my">平面坐标y（单位：meter）</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns></returns>
        public CoordinateInfo MetersToPixels(double mx, double my, int zoom)
        {
            CoordinateInfo retval = new CoordinateInfo();
            try
            {
                var res = Resolution(zoom);
                retval.X = (mx + this.originShift) / res;
                retval.Y = (my + this.originShift) / res;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 像素坐标转地图瓦片行列号
        /// </summary>
        /// <param name="px">像素坐标x</param>
        /// <param name="py">像素坐标y</param>
        /// <returns></returns>
        public TileIndexInfo PixelsToTile(double px, double py)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                retval.Column = (int)(Math.Ceiling(Convert.ToDouble(px / this.tileSize)) - 1);
                retval.Row = (int)(Math.Ceiling(Convert.ToDouble(py / this.tileSize)) - 1);
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// web墨卡托坐标转地图瓦片行列号
        /// </summary>
        /// <param name="mx">平面坐标x（单位：meter）</param>
        /// <param name="my">平面坐标y（单位：meter）</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns></returns>
        public TileIndexInfo MetersToTile(double mx, double my, int zoom)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                var p = this.MetersToPixels(mx, my, zoom);
                retval = this.PixelsToTile(p.X, p.Y);
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// WGS84地理坐标转转地图瓦片行列号（TMS服务）
        /// </summary>
        /// <param name="lat">纬度</param>
        /// <param name="lon">经度</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns>地图瓦片行列号</returns>
        public TileIndexInfo LatLonToTileTMS(double lat, double lon, int zoom)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                var m = this.LatLonToMeters(lat, lon);
                retval = this.MetersToTile(m.X, m.Y, zoom);
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// WGS84地理坐标转转地图瓦片行列号（XYZ服务）
        /// </summary>
        /// <param name="lat">纬度</param>
        /// <param name="lon">经度</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns>地图瓦片行列号</returns>
        public TileIndexInfo LatLonToTileXYZ(double lat, double lon, int zoom)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                var m = this.LatLonToMeters(lat, lon);
                retval = this.MetersToTile(m.X, m.Y, zoom);
                retval.Row = (int)Math.Pow(2, zoom) - retval.Row - 1;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 根据对应点行列号和地图缩放层级计算对应web墨卡托坐标系下四至坐标信息
        /// </summary>
        /// <param name="tx">column</param>
        /// <param name="ty">row</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns></returns>
        public GeoExtent TileBounds(int tx, int ty, int zoom)
        {
            GeoExtent retval = new GeoExtent();
            try
            {
                var min = this.PixelsToMeters(tx * this.tileSize, ty * this.tileSize, zoom);
                var max = this.PixelsToMeters((tx + 1) * this.tileSize, (ty + 1) * this.tileSize, zoom);
                retval = new GeoExtent() { North = max.Y, South = min.Y, East = max.X, West = min.X };
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 根据对应点行列号和地图缩放层级计算对应WGS84坐标系下四至坐标信息
        /// </summary>
        /// <param name="tx">column</param>
        /// <param name="ty">row</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns></returns>
        public GeoExtent TileLatLonBounds(int tx, int ty, int zoom)
        {
            GeoExtent retval = new GeoExtent();
            try
            {
                var bounds = this.TileBounds(tx, ty, zoom);
                var min = this.MetersToLatLon(bounds.West, bounds.South);
                var max = this.MetersToLatLon(bounds.East, bounds.North);
                retval = new GeoExtent() { North = max.Lat, South = min.Lat, East = max.Lng, West = min.Lng };
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// TMS协议行列号转Google地图行列号（TMS转XYZ）
        /// </summary>
        /// <param name="tx">column</param>
        /// <param name="ty">row</param>
        /// <param name="zoom">地图缩放层级zoom</param>
        /// <returns></returns>
        public TileIndexInfo TMSTileToGoogleTile(int tx, int ty, int zoom)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                retval.Column = tx;
                retval.Row = Convert.ToInt32((Math.Pow(2, zoom) - 1) - ty);
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 地图行列号转微软地图QuadTree（四叉树信息）
        /// </summary>
        /// <param name="tx"></param>
        /// <param name="ty"></param>
        /// <param name="zoom"></param>
        /// <returns></returns>
        public string QuadTree(int tx, int ty, int zoom)
        {
            string retval = "";
            try
            {

                ty = ((1 << zoom) - 1) - ty;
                for (var i = zoom; i >= 1; i--)
                {
                    var digit = 0;

                    var mask = 1 << (i - 1);

                    if ((tx & mask) != 0)
                        digit += 1;

                    if ((ty & mask) != 0)
                        digit += 2;

                    retval += digit;
                }

                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 根据微软地图QuadTree（四叉树信息）转行列号(TMS)
        /// </summary>
        /// <param name="quadtree"></param>
        /// <param name="zoom"></param>
        /// <returns></returns>
        public TileIndexInfo QuadTreeToTile(string quadtree, int zoom)
        {
            TileIndexInfo retval = new TileIndexInfo();
            try
            {
                var tx = 0;
                var ty = 0;

                for (var i = zoom; i >= 1; i--)
                {
                    var ch = quadtree[zoom - i];
                    var mask = 1 << (i - 1);

                    var digit = ch - '0';

                    if (Convert.ToBoolean(digit & 1))
                        tx += mask;

                    if (Convert.ToBoolean(digit & 2))
                        ty += mask;
                }

                ty = ((1 << zoom) - 1) - ty;
                retval.Column = tx;
                retval.Row = ty;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

        /// <summary>
        /// 经纬度转微软地图QuadTree（四叉树信息）
        /// </summary>
        /// <param name="lat">纬度</param>
        /// <param name="lon">经度</param>
        /// <param name="zoom">地图层级</param>
        /// <returns></returns>
        public string LatLonToQuadTree(double lat, double lon, int zoom)
        {
            string retval = "";
            try
            {

                var m = this.LatLonToMeters(lat, lon);
                var t = this.MetersToTile(m.X, m.Y, zoom);

                retval = this.QuadTree(Convert.ToInt32(t.Column), Convert.ToInt32(t.Row), zoom);

                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }
    }
}

• 真实测试

 

var lon = 108.35656;
var lat = 22.79632;
var zoom = 16;

var g = new WebMercator();
TileIndexInfo addressXYZ =g.LatLonToTileXYZ(lat, lon, zoom);
TileIndexInfo addressTMS = g.LatLonToTileTMS(lat, lon, zoom);
string quadTree = g.LatLonToQuadTree(lat, lon, zoom);
TileIndexInfo tempAddress =g.QuadTreeToTile(quadTree, zoom);

Console.WriteLine(tempAddress.ToString());
Console.WriteLine(addressXYZ.ToString());
Console.WriteLine(addressTMS.ToString());


Console.WriteLine("https://t1.ssl.ak.dynamic.tiles.virtualearth.net/comp/ch/"+ quadTree+"?mkt=en-US&it=A,G,L&src=o&og=1278&n=z");
//{IpServer}是geoserver所在服务器ip+port 
Console.WriteLine("http://{IpServer}/geoserver/gwc/service/wmts/rest/guangxi:img_c/raster/EPSG:900913/EPSG:900913:" +zoom + "/" + addressXYZ.Row + "/" + addressXYZ.Column + "?format=image/png");
Console.WriteLine("http://{IpServer}/geoserver/gwc/service/tms/1.0.0/guangxi%3Aimg_c@EPSG%3A900913@png/"
            + zoom + "/" + addressTMS.Column + "/" + addressTMS.Row + ".png");
Console.WriteLine("http://t2.tianditu.gov.cn/DataServer?T=img_w&x=" + addressXYZ.Column + "&y=" + addressXYZ.Row + "&l=" + zoom + "&tk=71d235688d12e7b0a230c5cc916b8d66");