wgs84座標系與gcj02座標系轉換誤差分佈圖 | Mapping the Error in Transformation between WGS84 and GCJ02 Coordinations

天靖居士發表於2017-10-12

國際上通用的是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中做的,如下圖所示,顏色越深代表誤差越大,越淺則誤差越小。
從圖中可以看出中國的三級階梯分佈。
wgs84座標系與gcj02座標系轉換誤差分佈圖 | Mapping the Error in Transformation between WGS84 and GCJ02 Coordinations

相關文章