C# 经纬度坐标的精度及WGS84(谷歌)、GCJ02(高德)、BD09(百度)坐标相互转换（含高精度转换）

1. 概述

WGS-84坐标系（World Geodetic System一1984 Coordinate System）是一种国际上采用的地心坐标系，GCJ-02是由中国国家测绘局（G表示Guojia国家，C表示Cehui测绘，J表示Ju局）制订的地理信息系统的坐标系统，BD09坐标系(Baidu Coordinate System)是百度地图使用的一种坐标系，它基于GCJ-02坐标系进行了加密偏移，适用于百度地图定位和导航服务。谷歌采用的是WGS-84坐标系，高德采用的是GCJ-02坐标系，也称火星坐标系。它们都是用经纬度表示其坐标。

2. 精度

我们通常在手机上见到的经纬度是用度、分、秒表示的，但是在以上三种坐标体系中常用的是以为度为单位的整数加6位小数的数值表示的。有时我们看到经纬度小数点后面很多位，似乎精度很高，这有意义吗？实际在我们生活中，GPS、北斗导航或定位时精确到几米就够用了，因为在几米范围内，你应该是能看到你要找的地方的。而“米”这个精度，正是经纬度小数点后的第5位，第6位则是分米，那么以上三种坐标用6位小数的精度，是能满足我们日常生活需要了。如果是要更高精度定位、室内定位等，单用GPS是满足不了的，需要使用UWB、差分等定位技术。
纬度小数点后面的每位精度是多少？一度的精度是多少？其实很好计算，我们知道，纬度是平行于赤道的圆，所以每度的长度是一样的。虽然地球是个不规则的椭圆球体，但我们可以假设地图是个标准球体，已知地球半径R是6371公里，因此垂直于赤道穿越地心圆周长是2πR=40009.88公里，除以360º，则1º=111.14公里，即纬度差1º，则差出约111公里，那么0.1º自然是11.1公里，0.01是1.11公里，0.001是111米，0.0001是11.1米，0.00001也就是1.11米了，从这儿可以看出小数点第5位其精度是米了，依次类推，小数点后第6位是分米，第7位是厘米，第8位是毫米，再往后去要求精度，实际意义已经不大了。
经度呢？道理一样，但不同的是经度构成的圆都是穿越地心的，赤道的一度的长度最长，南北极最短，不同的纬度上，经度一度的长度是不一样的，我们可以通过几何学去求解。假设纬度是θ，在此纬度的纬线构成的平行于赤道的圆半径则是地球半径乘以Cosθ，周长自然是2πRCosθ，再除以360º，即可得出此纬度上每度经度的长度，这个大家可以自己去根据自己所处的纬度实际去算一下。我们以赤道最大的长度来看，经度也是1º=111.14公里，跟纬度一样，其后面的小数点位代表的精度也跟纬度是一样的，就不再赘述了。
综上所述，经纬度的精度我们可以粗略地认为是一样的，即：
小数点后1位‌：约11114米（11公里）
‌小数点后2位‌：约1111米（1公里）
‌小数点后3位‌：约111米
‌小数点后4位‌：约11米
‌小数点后5位‌：约1米
‌小数点后6位‌：约0.1米（10厘米）
‌小数点后7位‌：约0.01米（1厘米）
‌小数点后8位‌：约0.001米（1毫米）

3. 经纬度坐标WGS84、GCJ-02、BD-09坐标相互转换代码

public class WGS84
{private const double pi = Math.PI;private const double xPi = Math.PI * 3000.0d / 180.0d;/***  判断经纬度是否超出中国境内*/public bool OutOfChina(double latitude, double longitude){if (longitude < 72.004 || longitude > 137.8347 || latitude < 0.8293 || latitude > 55.8271)return true;return false;}private double Offset_gcj02_latitude(double x, double y){double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.Sqrt(Math.Abs(x));ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0;ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0;return ret;}private double Offset_gcj02_longitude(double x, double y){double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.Sqrt(Math.Abs(x));ret += (20.0 * Math.Sin(6.0 * x * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0;ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0;ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0 * pi)) * 2.0 / 3.0;return ret;}private void Gcj02_delta(double lat, double lon, ref double dLat, ref double dLon){const double a = 6378245.0;                 // 长半轴const double ee = 0.00669342162296594323;   // 扁率double radLat = lat / 180.0 * pi;double magic = Math.Sin(radLat);magic = 1 - ee * magic * magic;double sqrtMagic = Math.Sqrt(magic);dLat = (Offset_gcj02_latitude(lon - 105.0, lat - 35.0) * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);dLon = (Offset_gcj02_longitude(lon - 105.0, lat - 35.0) * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);}/// <summary>/// WGS84坐标转换为GCJ02坐标/// </summary>/// <param name="wgsLat">WGS84纬度坐标(度)</param>/// <param name="wgsLon">WGS84经度坐标(度)</param>/// <param name="gcjLat">返回GCJ02纬度坐标(度)</param>/// <param name="gcjLon">返回GCJ02经度坐标(度)</param>public void Wgs84_to_gcj02(double wgsLat, double wgsLon, ref double gcjLat, ref double gcjLon){if (OutOfChina(wgsLat, wgsLon)){gcjLat = wgsLat;gcjLon = wgsLon;return;}double lat_offset = 0, lon_offset = 0;Gcj02_delta(wgsLat, wgsLon, ref lat_offset, ref lon_offset);gcjLat = wgsLat + lat_offset;gcjLon = wgsLon + lon_offset;}/// <summary>/// WGS84坐标转换为BD09坐标/// </summary>/// <param name="wgsLat">WGS84纬度坐标(度)</param>/// <param name="wgsLon">WGS84经度坐标(度)</param>/// <param name="bdLat">返回BD09纬度坐标</param>/// <param name="bdLon">返回BD09经度坐标</param>public void Wgs84_to_bd09(double wgsLat, double wgsLon, ref double bdLat, ref double bdLon){double gcjLat = 0.0d, gcjLon = 0.0d;Wgs84_to_gcj02(wgsLat, wgsLon, ref gcjLat, ref gcjLon);Gcj02_to_bd09(gcjLat, gcjLon, ref bdLat, ref bdLon);}/// <summary>/// GCJ02坐标转换为WGS84坐标(方法1：精度约为米，理论误差10米之内)/// </summary>/// <param name="gcjLat">GCJ02纬度坐标</param>/// <param name="gcjLon">GCJ02经度坐标</param>/// <param name="wgsLat">返回WGS84纬度坐标</param>/// <param name="wgsLon">返回WGS84经度坐标</param>public void Gcj02_to_wgs84(double gcjLat, double gcjLon, ref double wgsLat, ref double wgsLon){if (OutOfChina(gcjLat, gcjLon)){wgsLat = gcjLat;wgsLon = gcjLon;return;}const double a = 6378245.0;                 // 长半轴const double ee = 0.00669342162296594323;   // 扁率double dLat = Offset_gcj02_latitude(gcjLon - 105.0, gcjLat - 35.0);double dLng = Offset_gcj02_longitude(gcjLon - 105.0, gcjLat - 35.0);double radlat = gcjLat / 180.0 * pi;double magic = Math.Sin(radlat);magic = 1 - ee * magic * magic;double sqrtmagic = Math.Sqrt(magic);dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi);dLng = (dLng * 180.0) / (a / sqrtmagic * Math.Cos(radlat) * pi);double mgLat = gcjLat + dLat;double mgLng = gcjLon + dLng;wgsLon = gcjLon * 2 - mgLng;wgsLat = gcjLat * 2 - mgLat;}/// <summary>/// GCJ02坐标转换为WGS84坐标(方法2：高精度，几乎相等)/// </summary>/// <param name="gcjLat">GCJ02纬度坐标</param>/// <param name="gcjLon">GCJ02经度坐标</param>/// <param name="wgsLat">返回WGS84纬度坐标</param>/// <param name="wgsLon">返回WGS84经度坐标</param>/// <param name="threshold">精度阀值，默认为1e-17，转换后几乎相等</param>public void Gcj02_to_wgs84_2(double gcjLat, double gcjLon, ref double wgsLat, ref double wgsLon, double threshold = 1e-17){if (OutOfChina(gcjLat, gcjLon)){wgsLat = gcjLat;wgsLon = gcjLon;return;}const int maxIterations = 1000;double[] latBounds = { gcjLat - 0.5, gcjLat + 0.5 };double[] lonBounds = { gcjLon - 0.5, gcjLon + 0.5 };// 使用二分法逼近纬度double currentLat = (latBounds[0] + latBounds[1]) / 2;double currentLon = (lonBounds[0] + lonBounds[1]) / 2;int iteration = 0;double dLat, dLon;double tmpLat = 0.0d;double tmpLon = 0.0d;do{currentLat = (latBounds[0] + latBounds[1]) / 2;currentLon = (lonBounds[0] + lonBounds[1]) / 2;Wgs84_to_gcj02(currentLat, currentLon, ref tmpLat, ref tmpLon);dLat = tmpLat - gcjLat;dLon = tmpLon - gcjLon;if (dLat > 0){latBounds[1] = currentLat;}else{latBounds[0] = currentLat;}if (dLon > 0){lonBounds[1] = currentLon;}else{lonBounds[0] = currentLon;}iteration++;} while ((Math.Abs(dLat) > threshold || Math.Abs(dLon) > threshold) && iteration < maxIterations);wgsLat = currentLat;wgsLon = currentLon;}/// <summary>/// GCJ02坐标转换为BD09坐标/// </summary>/// <param name="gcjLat">GCJ02纬度坐标</param>/// <param name="gcjLon">GCJ02经度坐标</param>/// <param name="bdLat">返回BD09纬度坐标</param>/// <param name="bdLon">返回BD09经度坐标</param>public void Gcj02_to_bd09(double gcjLat, double gcjLon, ref double bdLat, ref double bdLon){double z = Math.Sqrt(gcjLon * gcjLon + gcjLat * gcjLat) + 0.00002 * Math.Sin(gcjLat * xPi);double theta = Math.Atan2(gcjLat, gcjLon) + 0.000003 * Math.Cos(gcjLon * xPi);bdLon = z * Math.Cos(theta) + 0.0065;bdLat = z * Math.Sin(theta) + 0.006;}/// <summary>/// BD09坐标转换为GCJ02坐标(方法1，精度约为分米，理论误差小于1米)/// </summary>/// <param name="bdLat">BD09纬度坐标</param>/// <param name="bdLon">BD09经度坐标</param>/// <param name="gcjLat">返回GCJ02纬度坐标</param>/// <param name="gcjLon">返回GCJ02经度坐标</param>public void Bd09_to_gcj02(double bdLat, double bdLon, ref double gcjLat, ref double gcjLon){double x = bdLon - 0.0065;double y = bdLat - 0.006;double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * xPi);double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * xPi);gcjLon = z * Math.Cos(theta);gcjLat = z * Math.Sin(theta);}/// <summary>/// BD09坐标转换为GCJ02坐标(方法2，高精度，转换后几乎相等)/// </summary>/// <param name="bdLat">BD09纬度坐标</param>/// <param name="bdLon">BD09经度坐标</param>/// <param name="gcjLat">返回GCJ02纬度坐标</param>/// <param name="gcjLon">返回GCJ02经度坐标</param>/// <param name="threshold">精度阀值，默认为1e-17，转换后几乎相等</param>public void Bd09_to_gcj02_2(double bdLat, double bdLon, ref double gcjLat, ref double gcjLon, double threshold = 1e-17){double lon = bdLon;double lat = bdLat;double x = lon - 0.0065;double y = lat - 0.006;gcjLon = x;gcjLat = y;double tempLat = 0.0d;double tempLon = 0.0d;Gcj02_to_bd09(gcjLat, gcjLon, ref tempLat, ref tempLon);double dx = tempLon - lon;double dy = tempLat - lat;while (Math.Abs(dx) > threshold || Math.Abs(dy) > threshold)//回归{gcjLon -= dx;gcjLat -= dy;Gcj02_to_bd09(gcjLat, gcjLon, ref tempLat, ref tempLon);dx = tempLon - lon;dy = tempLat - lat;}}/// <summary>/// BD09坐标转换为WGS84坐标(方法1，精度约为米或十几米左右)/// </summary>/// <param name="bdLat">BD09纬度坐标</param>/// <param name="bdLon">BD09经度坐标</param>/// <param name="wgsLat">返回WGS84纬度坐标</param>/// <param name="wgsLon">返回WGS84经度坐标</param>public void Bd09_to_wgs84(double bdLat, double bdLon, ref double wgsLat, ref double wgsLon){double gcjLat = 0.0d, gcjLon = 0.0d;Bd09_to_gcj02(bdLat, bdLon, ref gcjLat, ref gcjLon);Gcj02_to_wgs84(gcjLat, gcjLon, ref wgsLat, ref wgsLon);}/// <summary>/// BD09坐标转换为WGS84坐标(方法2，高精度，转换后几乎相等)/// </summary>/// <param name="bdLat">BD09纬度坐标</param>/// <param name="bdLon">BD09经度坐标</param>/// <param name="wgsLat">返回WGS84纬度坐标</param>/// <param name="wgsLon">返回WGS84经度坐标</param>/// <param name="gcj_threshold">转换GCJ02精度阀值，默认为1e-17，转换后几乎相等</param>/// <param name="wgs_threshold">转换WGS84精度阀值，默认为1e-17，转换后几乎相等</param>public void Bd09_to_wgs84_2(double bdLat, double bdLon, ref double wgsLat, ref double wgsLon, double gcj_threshold = 1e-17, double wgs_threshold = 1e-17){double gcjLat = 0.0d, gcjLon = 0.0d;Bd09_to_gcj02_2(bdLat, bdLon, ref gcjLat, ref gcjLon, gcj_threshold);Gcj02_to_wgs84_2(gcjLat, gcjLon, ref wgsLat, ref wgsLon, wgs_threshold);}/// <summary>/// WGS84大地坐标转换为ECEF空间直角坐标/// </summary>/// <param name="lat">WGS84纬度坐标</param>/// <param name="lon">WGS84经度坐标</param>/// <param name="alt">WGS84高度坐标</param>/// <param name="x">返回ECEF空间X坐标</param>/// <param name="y">返回ECEF空间Y坐标</param>/// <param name="z">返回ECEF空间Z坐标</param>public void Wgs84_to_ecef(double lat, double lon, double alt, ref double x, ref double y, ref double z){const double a = 6378137.0;             // WGS-84 椭球体的长半轴（赤道半径）const double f = 1 / 298.257223563;     // 地球扁率const double e2 = 2 * f - f * f;        // 0.00669437999014;  WGS-84 椭球体的第一偏心率的平方double radLat = lat * pi / 180.0;     // 将纬度转换为弧度double radLon = lon * pi / 180.0;     // 将经度转换为弧度double N = a / Math.Sqrt(1 - e2 * Math.Sin(radLat) * Math.Sin(radLat));  // 椭球曲率半径x = (N + alt) * Math.Cos(radLat) * Math.Cos(radLon);y = (N + alt) * Math.Cos(radLat) * Math.Sin(radLon);z = (N * (1 - e2) + alt) * Math.Sin(radLat);}/// <summary>/// ECEF空间直角坐标转换为WGS84大地坐标/// </summary>/// <param name="x">ECEF空间X坐标</param>/// <param name="y">ECEF空间Y坐标</param>/// <param name="z">ECEF空间Z坐标</param>/// <param name="lat">返回WGS84纬度坐标</param>/// <param name="lon">返回WGS84经度坐标</param>/// <param name="alt">返回WGS84高度坐标</param>public void Ecef_to_wgs84(double x, double y, double z, ref double lat, ref double lon, ref double alt){const double a = 6378137.0;                 // WGS-84 椭球体的长半轴（赤道半径）const double f = 1.0 / 298.257223563;       // 扁率const double e2 = 2 * f - f * f;            // 0.00669437999014;  WGS-84 椭球体的第一偏心率的平方const double b = a * (1 - f);               // 6356752.314245; WGS-84 椭球体的短半轴（极半径）double p = Math.Sqrt(x * x + y * y);lon = Math.Atan2(y, x);    // 计算经度// 迭代计算纬度double theta = Math.Atan2(z * a, p * b);double sinTheta = Math.Sin(theta);double cosTheta = Math.Cos(theta);lat = Math.Atan2(z + e2 * b * sinTheta * sinTheta * sinTheta,p - e2 * a * cosTheta * cosTheta * cosTheta);double N = a / Math.Sqrt(1 - e2 * Math.Sin(lat) * Math.Sin(lat));alt = p / Math.Cos(lat) - N;lat = lat * 180.0 / pi;   // 转换为度数lon = lon * 180.0 / pi;   // 转换为度数}}

4. 调用

因为返回经纬度坐标的变量使用的是ref关键字，所以必须先声明纬度、经度变量，例如WGS84坐标转GCJ-02坐标：

private WGS84 WGS = new WGS84();double WGSLat = 37.786787025540519d;
double WGSLng = 112.54513755035273d;double GCJLat = 0.0d;
double GCJLng = 0.0d;WGS.Wgs84_to_gcj02(WGSLat, WGSLng, ref GCJLat, ref GCJLng);string GCJLatString = GCJLat.ToString("R");
string GCJLng = GCJLng.ToString("R");

注意：代码中ToString()方法中用的是“R”参数，这样能把小数点后的位数全部显示出来，如果不加参数，只能显示小数点后13位。

其他坐标转换方法大致相同，参考上面的例子即可。

5. 说明

（1）转换代码中，转换方向是WGS84转GCJ-02，GCJ-02再转BD-09，回归转换时是BD-09转GCJ-02，GCJ-02再转WGS84。
（2）回归转换增加了高精度转换方法：Gcj02_to_wgs84_2()以及Bd09_to_gcj02_2()，由这两个方法组合成了第三个高精度转换方法Bd09_to_wgs84_2()。
(3) 回归的非高精度的转换方法足能满足我们正常应用，我测试过几次，BD-09转GCJ-02误差在小数点五六位，转WGS84一般在小数点第五位，也有在第四位上，但也就误差在1左右，即十几米。
（4）回归高精度转换方法，误差在小数点后15位左右，转换精度几乎一样。代码中的方法采用普通回归方法还是高精度回归方法，大家根据自己需求决定。
(5) GCJ-02转BD-09以及BD-09转GCJ-02（非高精度）使用的是xPi，而不是Pi（π），xPi = Pi * 3000.0 / 180.0，这跟其他方法中Pi是有区别的。