中国地图坐标(GCJ-02)偏移算法破解小史

2006年，Google开始与AutoNavi合作使用后者所提供的中国地图。这应该是外企首次接触到这个问题。

从2009年开始，中国地图的坐标偏移开始为外界所知。Garmin的用户发现在美国购买的GPS到了中国几乎无法使用，而在中国购买的Garmin产品则没有问题。Google Maps API的使用者发现兴趣点无法被准确标注在中国地图上。更有意思的是，有用户反复就此报告bug给Google，却从未得到任何回应。类似的，Garmin也声称自己没有解决方案，建议客户在需要的情况下在中国境内购买GPS设备。

于此同时，各路豪杰开始尝试破解这种偏移算法。其中有两条路径值得注意：

2010年1月，网友wuyongzheng发现：

I accidentally found the Chinese version of Google Map ditu.google.com to be able to correlate satellite image with map, and it gives the amount of deviation for any location in China. This URL queries the deviation of 34.29273N,108.94695E (Xi’an): http://ditu.google.com/maps/vp?spn=0.001,0.001&t=h&z=18&vp=$34.29273,108.94695 (seems it’ doesn’t work now)

有了足够的数据，wuyongzheng建议使用回归算法来逼近这个偏移算法：https://wuyongzheng.wordpress.com/2010/01/22/china-map-deviation-as-a-regression-problem/

在此之前的尝试都是零星的，针对个别城市的。wuongzheng的这个建议算是在全面系统地解决这个问题上迈出了第一步。

2013年5月，Maxime Guilbot根据这个建议得到4-5米精度的逼近：

https://github.com/maxime/ChinaMapDeviation

2013年10月，wuyongzheng自己进行了回归，得到如下结果：

http://wuyongzheng.github.io/china-map-deviation/paper.html

Maxime Guibot和wuyongzheng的回归结果基本代表了在黑暗中摸索的最佳结果，因此得到了广泛的注意和应用。

在另一条路径上，2010年4月，emq project增加了一个文件，Converter.java：

http://emq.googlecode.com/svn/emq/src/Algorithm/Coords/Converter.java

这段代码可以以很高的精度把WGS-84坐标转换到GCJ-02坐标。

2013年2月，这段代码被网友coolypf注意到，整理后用到了他自己的项目中：

https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#353936

其中的关键代码值得贴在这里：

        const double pi = 3.14159265358979324;

        //
        // Krasovsky 1940
        //
        // a = 6378245.0, 1/f = 298.3
        // b = a * (1 - f)
        // ee = (a^2 - b^2) / a^2;
        const double a = 6378245.0;
        const double ee = 0.00669342162296594323;

        //
        // World Geodetic System ==> Mars Geodetic System
        public static void transform(double wgLat, double wgLon, out double mgLat, out double mgLon)
        {
            if (outOfChina(wgLat, wgLon))
            {
                mgLat = wgLat;
                mgLon = wgLon;
                return;
            }
            double dLat = transformLat(wgLon - 105.0, wgLat - 35.0);
            double dLon = transformLon(wgLon - 105.0, wgLat - 35.0);
            double radLat = wgLat / 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);
            dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi);
            mgLat = wgLat + dLat;
            mgLon = wgLon + dLon;
        }

2013年3月，coolypf在自己的博客中介绍了这一段代码：

http://blog.csdn.net/coolypf/article/details/8686588

2014年9月，wuyongzheng注意到了coolypf的项目。至此，两条路径合流，坐标偏移问题基本得到了完美解决。

从上面的代码可以看出，相对于WGS－84，GCJ－02一方面采用了不同的参考椭球体(SK-42, Krasovsky。应该属于前苏联影响的遗留)，另一方面引入了高频非线性偏移。