在做项目的时候，使用OpenLayers地图来绘制范围和工程的关系，出现了一个很有意思的情况。通常来说计算坐标点是否包含在圆内是通过计算两点之间的距离和所画圆的半径来进行比较，基于这种计算，圆外却出现了乱入工程，如下图

疑问来了，到底是计算错了呢还是范围错了呢，通过代码计算反复验证后，结果都是对的，圆外标红的位置通过Openlayers内置的计算结果也是对的，为什么没有包含在圆里呢？

这个时候，我们先来说一下墨卡托投影，问题就会迎刃而解。

什么是墨卡托投影？
墨卡托（Mercator）投影，又名“等角正轴圆柱投影”，荷兰地图学家墨卡托（Mercator）在1569年拟定，假设地球被围在一个中空的圆柱里，其赤道与圆柱相接触，然后再假想地球中心有一盏灯，把球面上的图形投影到圆柱体上，再把圆柱体展开，这就是一幅标准纬线为零度（即赤道）的“墨卡托投影”绘制出的世界地图 。从球到平面，有转换公式，感兴趣的可以去查看（地址如下，公式和过程很全乎） 

 https://blog.csdn.net/cleverysm/article/details/2121633

为什么选择墨卡托投影？
      墨卡托投影的“等角”特性，说白了就是能够保证对对象外观的的正确性和属性（方向，位置）的正确性，在使用时候不会出错。另外墨卡托投影的“圆柱”特性，保证了纬线经线都是平行直线，并且相互垂直。而且经线间隔是相同的，纬线间隔从赤道向两级逐渐增大。但是，“等角”不可避免的带来的面积的变形，尤其是两极地区，格陵兰岛比实际面积扩大了N倍（后知后觉，以前一直觉得它面积好大），南极洲也是如此。可以对比一下，第一张为经纬度投影，纬度越高，形变越明显，第二张为墨卡托投影 ， 纵向距离也是随着纬度增大而变长， 横向变大，同时纵向也变大，而且变化比例接近（因为0.0000X纬度引起可忽略），结果就是只把一个图形“原样放大”了,而形状却没有变化。 值得注意的是，web墨卡托虽然形状没变，但是高纬度地区的面积比真实同样放大了很多倍，在中国地图上选取最南端的三亚和最北端的漠河，虽然形状没有改变，却在地图上的面积相差好几倍，“大漠河，小三亚”。

    图1 经纬度投影

  图2 墨卡托投影

Web墨卡托投影坐标系：
以整个世界范围，赤道作为标准纬线，本初子午线作为中央经线，两者交点为坐标原点，向东向北为正，向西向南为负。
X轴：由于赤道半径为6378137米，则赤道周长为2*PI*r = 2*20037508.3427892，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。
Y轴：由墨卡托投影的公式可知，当纬度φ接近两极，即90°时，y值趋向于无穷。Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间，做成正方形。
       因此在投影坐标系（米）下的范围是：最小(-20037508.3427892, -20037508.3427892 )到最大 (20037508.3427892, 20037508.3427892)。
对应的地理坐标系:
      简单来说，投影坐标系(PROJCS)是平面坐标系，以米为单位；而地理坐标系(GEOGCS)是椭球面坐标系，以经纬度为单位。
经度：可取全球范围：[-180,180]。
纬度：我们都知道，纬度不可能到达90°，墨卡托投影做成的正方形取值-20037508.3427892，经过反计算，可得到纬度85.05112877980659。因此纬度取值范围是[-85.05112877980659，85.05112877980659]。 因此，地理坐标系（经纬度）对应的范围是：最小(-180,-85.05112877980659)，最大(180, 85.05112877980659)。
       除Google地图外，Bing地图、MapABC地图、百度地图、天地图等在线地图服务均采用了Web墨卡托投影，ESRI的Online地图也有使用此地图投影。除了在切片时候减少图片数量以外，也为了方便使用，均采用了Web 墨卡托投影方式，但在地图显示上面为了易于观看，还是转换成了经纬度的方式显示。

      在WebGIS中，经常有已知WGS84(EPSG:4326)经纬度[lng, lat]，在Canvas2D上要显示该点就要换算成墨卡托投影(EPSG:3857)，其实只要知道地球变成平面地图时的转换比率，那就能通过经纬度求出墨卡托的值。Openlayers以及Mapbox/sphericalmercator都有转换函数。

     回到开始的问题，展示出来的圆形是否是正确的呢，其实实际上应该是展现为椭圆形（本质上还是画圆）

附上圆形计算距离方法和多变形计算距离方法：

 public async Task<double> GetRealDistance(double sLng, double sLat, double eLng, double eLat)
        {
              double distance = Math.Round(
               6370.856 * 2 * Math.Asin(
                    Math.Sqrt(
                        Math.Pow(Math.Sin((sLat * Math.PI / 180 - eLat * Math.PI / 180)) / 2, 2) +
                        Math.Cos(sLat * Math.PI / 180) * Math.Cos(eLat * Math.PI / 180) *
                        Math.Pow(Math.Sin((sLng * Math.PI / 180 - eLng * Math.PI / 180)) / 2, 2)
                    )
                ) * 1000
            );          
            return distance;
        }

        /// <summary>
        /// 判断点是否在多边形内或多边形上（光线追踪算法简单应用）
        /// </summary>
        /// <param name="ALon">经度</param>
        /// <param name="ALat">纬度</param>
        /// <param name="Points">多边形边界点集合</param>
        /// <returns></returns>
        public static bool IsPtInPoly(double ALon, double ALat, Point[] Points)
        {
            int iSum, iCount, iIndex;
            double dLon1 = 0, dLon2 = 0, dLat1 = 0, dLat2 = 0, dLon;
            if (Points.Length < 3)
            {
                return false;
            }
            iSum = 0;
            iCount = Points.Length;
            for (iIndex = 0; iIndex < iCount; iIndex++)
            {
                if (ALon == Points[iIndex].getX() && ALat == Points[iIndex].getY())  //A点在多边形上    
                    return true;

                if (iIndex == iCount - 1)
                {
                    dLon1 = Points[iIndex].getX();
                    dLat1 = Points[iIndex].getY();
                    dLon2 = Points[0].getX();
                    dLat2 = Points[0].getY();
                }
                else
                {
                    dLon1 = Points[iIndex].getX();
                    dLat1 = Points[iIndex].getY();
                    dLon2 = Points[iIndex + 1].getX();
                    dLat2 = Points[iIndex + 1].getY();
                }

                //以下语句判断A点是否在边的两端点的纬度之间，在则可能有交点
                if (((ALat > dLat1) && (ALat < dLat2)) || ((ALat > dLat2) && (ALat < dLat1)))
                {
                    if (Math.Abs(dLat1 - dLat2) > 0)
                    {
                        //获取A点向左射线与边的交点的x坐标：
                        dLon = dLon1 - ((dLon1 - dLon2) * (dLat1 - ALat)) / (dLat1 - dLat2);
                        //如果交点在A点左侧，则射线与边的全部交点数加一：
                        if (dLon < ALon)
                        {
                            iSum++;
                        }
                        //如果相等，则说明A点在边上
                        if (dLon == ALon)
                            return true;
                    }
                }
            }
            if ((iSum % 2) != 0)
            {
                return true;
            }
            return false;
        }