wgs84坐标系与gcj02坐标系转换误差分布图 | Mapping the Error in Transformation between WGS84 and GCJ02 Coordinations

国际上通用的是wgs84坐标系，而我国对于境内的坐标进行了加密，采用了gcj02坐标系，或者称为火星坐标系。亢孟军老师带的一门课《多媒体电子地图设计》要求我们从wgs84坐标系转换为gcj02坐标系，再反算出wgs84坐标系并进行前后wgs84坐标系误差分析。在这里简单介绍一下方法。

生成等间距矩阵点

首先简单介绍一段代码，使用python写的，作用是生成等间距矩阵点。
下面是源代码的类 LocationDivide，可以直接将这段代码拷贝进行使用。

# Research region
class LocationDivide(object):
    def __init__(self, bound, size):
        # minLat,minLon,maxLat,maxLon
        self.minLat = float(str(bound).split(',')[0])
        self.minLon = float(str(bound).split(',')[1])
        self.maxLat = float(str(bound).split(',')[2])
        self.maxLon = float(str(bound).split(',')[3])
        self.size = size

    # Seperate block into blocks
    def compute_block(self):
        if (self.maxLat - self.minLat) % self.size == 0:
            lat_count = (self.maxLat - self.minLat) / self.size
        else:
            lat_count = (self.maxLat - self.minLat) / self.size + 1
        if (self.maxLon - self.minLon) % self.size == 0:
            lon_count = (self.maxLon - self.minLon) / self.size
        else:
            lon_count = (self.maxLon - self.minLon) / self.size + 1
        self.bounds = []
        lat_count = int(lat_count)
        lon_count = int(lon_count)

        try:
            for i in range(0, lat_count):
                for j in range(0, lon_count):
                    # maxLat,minLon,minLat,maxLon
                    minLat = self.minLat + i * self.size
                    minLon = self.minLon + j * self.size
                    maxLat = self.minLat + (i + 1) * self.size
                    if maxLat > self.maxLat:
                        maxLat = self.maxLat
                    maxLon = self.minLon + (j + 1) * self.size
                    if maxLon > self.maxLon:
                        maxLon = self.maxLon
                    # minLat,minLon,maxLat,maxLon
                    # set decimal
                    digit_number = 5
                    minLat = round(minLat, digit_number)
                    minLon = round(minLon, digit_number)
                    maxLat = round(maxLat, digit_number)
                    maxLon = round(maxLon, digit_number)
                    bound = "{0},{1},{2},{3}".format(minLat, minLon, maxLat, maxLon)
                    self.bounds.append(bound)
        except Exception as e:
            with open("e:log.txt", 'a') as log:
                log.writelines(e)
        return self.bounds

具体的使用方法非常简单，如下：

    # Set region bound and interval
    # minLat,minLon,maxLat,maxLon,interval

    region = "17,73,53,135"
    location = LocationDivide(region, 0.5)

    # Seperate grid into blocks

    location.compute_block()

即只要输入左下和右上角的经纬度坐标，再输入点之间的间隔，调用LocationDivide.compute_block()方法即可在划定的区域内，均匀生成等间距的采样点。需要使用时只需要调用LocationDivide.bounds变量即可，可以获得所有的点坐标，从左到右，从下到上依次排列。
这段代码非常实用，可以应用于多种场景，特别是需要生成格网时和需要均匀空间采样时可以使用。

wgs84坐标系转换为gcj02坐标系

在这里我们可以使用高德地图API进行数据的转换。具体可以参看其API
该API可以将wgs84坐标系、百度坐标系等坐标系的点转换为gcj02坐标系中点的值。
特别注意的是，坐标点经纬度小数不超过6位，且一次最多传入40对坐标点。

wgs84坐标系与gcj02坐标系之间的相互转换

虽然我们可以使用高德地图API将wgs84数据转换为gcj02坐标系的数据，但如何从gcj02反解出wgs84坐标是一个问题。当然网上有很多相应的资料，这里仅仅简单列举两条我的参考文献：
wgs84和gcj02 相互转换JAVA代码，包括我的代码也主要使用了这位博主的代码：
http://www.cnblogs.com/94cool/p/4266907.html
windows phone上wgs84转换成gcj02的C#代码，可以通过这个反推出gcj02计算wgs84的方法。
https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#353936
在这里还是附上我的代码，主要是使用了94cool博主的代码，修改为C#之后的代码，包括了两个类 Gps和PositionUtil：

public class Gps
        {
            double latitude { set; get; }
            double longitude { set; get; }
            public Gps(double latitude, double lontitude)
            {
                this.latitude = latitude;
                this.longitude = lontitude;
            }
            public Gps(string gps)
            {
                this.latitude = Convert.ToDouble(gps.Split(',')[1]);
                this.longitude = Convert.ToDouble(gps.Split(',')[0]);
            }
            public double getWgLat()
            {
                return this.latitude;
            }
            public double getWgLon()
            {
                return this.longitude;
            }
        }

public class PositionUtil
        {

            public static String BAIDU_LBS_TYPE = "bd09ll";

            public static double pi = 3.1415926535897932384626;
            public static double a = 6378245.0;
            public static double ee = 0.00669342162296594323;

            /**
             * 84 to 火星坐标系 (GCJ-02) World Geodetic System ==> Mars Geodetic System
             * 
             * @param lat
             * @param lon
             * @return
             */
            public static Gps gps84_To_Gcj02(double lat, double lon)
            {
                if (outOfChina(lat, lon))
                {
                    return null;
                }
                double dLat = transformLat(lon - 105.0, lat - 35.0);
                double dLon = transformLon(lon - 105.0, lat - 35.0);
                double radLat = lat / 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);
                double mgLat = lat + dLat;
                double mgLon = lon + dLon;
                return new Gps(mgLat, mgLon);
            }

            /**
             * * 火星坐标系 (GCJ-02) to 84 * * @param lon * @param lat * @return
             * */
            public static Gps gcj_To_Gps84(double lat, double lon)
            {
                Gps gps = transform(lat, lon);
                double lontitude = lon * 2 - gps.getWgLon();
                double latitude = lat * 2 - gps.getWgLat();
                return new Gps(latitude, lontitude);
            }

            /**
             * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 将 GCJ-02 坐标转换成 BD-09 坐标
             * 
             * @param gg_lat
             * @param gg_lon
             */
            public static Gps gcj02_To_Bd09(double gg_lat, double gg_lon)
            {
                double x = gg_lon, y = gg_lat;
                double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * pi);
                double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * pi);
                double bd_lon = z * Math.Cos(theta) + 0.0065;
                double bd_lat = z * Math.Sin(theta) + 0.006;
                return new Gps(bd_lat, bd_lon);
            }

            /**
             * * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 * * 将 BD-09 坐标转换成GCJ-02 坐标 * * @param
             * bd_lat * @param bd_lon * @return
             */
            public static Gps bd09_To_Gcj02(double bd_lat, double bd_lon)
            {
                double x = bd_lon - 0.0065, y = bd_lat - 0.006;
                double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * pi);
                double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * pi);
                double gg_lon = z * Math.Cos(theta);
                double gg_lat = z * Math.Sin(theta);
                return new Gps(gg_lat, gg_lon);
            }

            /**
             * (BD-09)-->84
             * @param bd_lat
             * @param bd_lon
             * @return
             */
            public static Gps bd09_To_Gps84(double bd_lat, double bd_lon)
            {

                Gps gcj02 = PositionUtil.bd09_To_Gcj02(bd_lat, bd_lon);
                Gps map84 = PositionUtil.gcj_To_Gps84(gcj02.getWgLat(),
                        gcj02.getWgLon());
                return map84;

            }

            public static Boolean outOfChina(double lat, double lon)
            {
                if (lon < 72.004 || lon > 137.8347)
                    return true;
                if (lat < 0.8293 || lat > 55.8271)
                    return true;
                return false;
            }

            public static Gps transform(double lat, double lon)
            {
                if (outOfChina(lat, lon))
                {
                    return new Gps(lat, lon);
                }
                double dLat = transformLat(lon - 105.0, lat - 35.0);
                double dLon = transformLon(lon - 105.0, lat - 35.0);
                double radLat = lat / 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);
                double mgLat = lat + dLat;
                double mgLon = lon + dLon;
                return new Gps(mgLat, mgLon);
            }

            public static double transformLat(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;
            }

            public static double transformLon(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;
            }

        }

调用方法直接使用PositionUtil.gcj_To_Gps84()等方法，传入相应的参数即可。

误差分布图

最后的误差分布图是在ArcGIS中做的，如下图所示，颜色越深代表误差越大，越浅则误差越小。
从图中可以看出中国的三级阶梯分布。