1. 程式人生 > >iOS開發中的火星座標系及各種座標系轉換演算法

iOS開發中的火星座標系及各種座標系轉換演算法

原文地址:http://m.oschina.net/blog/619183?ref=myread

其原理是這樣的:保密局開發了一個系統,能將實際的座標轉換成虛擬的座標。所有在中國銷售的數字地圖必須使用這個系統進行座標轉換之後方可上市。這是生產環節,這種電子地圖被稱為火星地圖。在使用環節,GPS終端裝置必須整合保密局提供的加密演算法(整合工作由保密局完成),把從GPS衛星那裡得到的座標轉換成虛擬座標,然後再去火星地圖上查詢,這樣就在火星座標系上完成了地圖的匹配。 所以大家所用的百度,高德等地圖定位準是偏差幾百米

名詞總結:

地球座標:指WGS84座標系統

火星座標:指使用國家保密外掛人為偏移後的座標


地球地圖:指與地球座標對應的客觀真實的地圖
火星地圖:指經過加密偏移後的,與火星座標對應的地圖

座標系轉換演算法

1.GCJ-02(火星座標系)和BD-09轉換

// GCJ-02 座標轉換成 BD-09 座標
+ (CLLocationCoordinate2D)MarsGS2BaiduGS:(CLLocationCoordinate2D)coordinate
{
    double x_pi = PI * 3000.0 / 180.0;
    double x = coordinate.longitude, y = coordinate.latitude;
    double z = sqrt(x * x + y * y) + 0.00002 * sin(y * x_pi);
    double theta = atan2(y, x) + 0.000003 * cos(x * x_pi);
    double bd_lon = z * cos(theta) + 0.0065;
    double bd_lat = z * sin(theta) + 0.006;
    return CLLocationCoordinate2DMake(bd_lat, bd_lon);
}

// BD-09 座標轉換成 GCJ-02 座標
+ (CLLocationCoordinate2D)BaiduGS2MarsGS:(CLLocationCoordinate2D)coordinate
{
    double x_pi = PI * 3000.0 / 180.0;
    double x = coordinate.longitude - 0.0065, y = coordinate.latitude - 0.006;
    double z = sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi);
    double theta = atan2(y, x) - 0.000003 * cos(x * x_pi);
    double gg_lon = z * cos(theta);
    double gg_lat = z * sin(theta);
    return CLLocationCoordinate2DMake(gg_lat, gg_lon);
}

2WGS-84(地球座標系)和BD-09(百度座標)轉換

// WGS-84 座標轉換成 BD-09 座標
+ (CLLocationCoordinate2D)WorldGS2BaiduGS:(CLLocationCoordinate2D)coordinate
{
    CLLocationCoordinate2D mars = [ALDGeocoder WorldGS2MarsGS:coordinate];
    CLLocationCoordinate2D baidu = [ALDGeocoder MarsGS2BaiduGS:mars];
    return baidu;
}

// BD-09 座標轉換成 WGS-84 座標
+ (CLLocationCoordinate2D)BaiduGS2WorldGS:(CLLocationCoordinate2D)coordinate
{
    CLLocationCoordinate2D mars = [ALDGeocoder BaiduGS2MarsGS:coordinate];
    CLLocationCoordinate2D world = [ALDGeocoder MarsGS2WorldGS:mars];
    return world;
}

3.WGS-84和sogou座標轉換

// WGS-84 座標轉換成 Sogou 座標
+ (CLLocationCoordinate2D)WorldGS2SogouGS:(CLLocationCoordinate2D)coordinate
{
    const double ee = 0.082271854224939184;
    double lon = coordinate.longitude;
    double lat = coordinate.latitude;
    double dlon = [ALDGeocoder rad:CLIP(lon, -360, 360)];
    double dlat = [ALDGeocoder rad:CLIP(lat, -90, 90)];
    dlon = 6378206.4 * dlon;
    double sinphi = sin(dlat);
    double temp1, temp2;
    if((temp1 = 1.0 + sinphi) == 0.0){
        dlat = -1000000000;
    }else if((temp2 = 1.0 - sinphi) == 0.0){
        dlat = 1000000000;
    }else{
        double esinphi = ee * sinphi;
        dlat = 3189103.2000000002 * log((temp1 / temp2) * pow((1.0 - esinphi) / (1.0 + esinphi), ee));
    }
    return CLLocationCoordinate2DMake(dlat, dlon);
}

// Sogou 座標轉換成 WGS-84 座標
+ (CLLocationCoordinate2D)SogouGS2WorldGS:(CLLocationCoordinate2D)coordinate
{
    const double ee = 1.5707963267948966;
    const double aa = 0.0033938814110493522;
    double lon = coordinate.longitude;
    double lat = coordinate.latitude;
    double dlon = lon / 6378206.4;
    double temp = -lat / 6378206.4;
    double chi;
    if(temp < -307){
        chi = ee;
    }else if(temp > 308){
        chi = -ee;
    }else{
        chi = ee - 2 * atan(exp(temp));
    }
    double chi2 = 2 * chi;
    double coschi2 = cos(chi2);
    double dlat = chi + sin(chi2) * (aa + coschi2 * (1.3437644537757259E-005 + coschi2 * (7.2964865099246009E-008 + coschi2 * 4.4551470401894685E-010)));
    double rlon = CLIP([ALDGeocoder deg:dlon], -360, 360);
    double rlat = CLIP([ALDGeocoder deg:dlat], -90, 90);
    return CLLocationCoordinate2DMake(rlat, rlon);
}

4火星座標和地球座標轉換

// World Geodetic System ==> Mars Geodetic System
+ (CLLocationCoordinate2D)WorldGS2MarsGS:(CLLocationCoordinate2D)coordinate
{
    // 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;
    
    if (outOfChina(coordinate.latitude, coordinate.longitude))
    {
        return coordinate;
    }
    double wgLat = coordinate.latitude;
    double wgLon = coordinate.longitude;
    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 = sin(radLat);
    magic = 1 - ee * magic * magic;
    double sqrtMagic = sqrt(magic);
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI);
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * PI);
    
    return CLLocationCoordinate2DMake(wgLat + dLat, wgLon + dLon);
}

// Mars Geodetic System ==> World Geodetic System
+ (CLLocationCoordinate2D)MarsGS2WorldGS:(CLLocationCoordinate2D)coordinate
{
    double gLat = coordinate.latitude;
    double gLon = coordinate.longitude;
    CLLocationCoordinate2D marsCoor = [ALDGeocoder WorldGS2MarsGS:coordinate];
    double dLat = marsCoor.latitude - gLat;
    double dLon = marsCoor.longitude - gLon;
    return CLLocationCoordinate2DMake(gLat - dLat, gLon - dLon);
}

5WGS-84 墨卡託座標轉換

//WGS-84 座標轉換成 墨卡託 座標
+ (CLLocationCoordinate2D)WorldGS2Mercator:(CLLocationCoordinate2D)coordinate
{
    double lon = coordinate.longitude*20037508.34/180;
    double lat = log(tan((90+coordinate.latitude)*M_PI/360))/(M_PI/180);
    lat = lat*20037508.34/180;
    return CLLocationCoordinate2DMake(lat, lon);
}

//墨卡託 座標轉換成 WGS-84 座標
+ (CLLocationCoordinate2D)Mercator2WorldGS:(CLLocationCoordinate2D)mercator
{
    double lon = mercator.longitude/20037508.34*180;
    double lat = mercator.latitude/20037508.34*180;
    lat = 180/M_PI*(2*atan(exp(lat*M_PI/180))-M_PI/2);
    return CLLocationCoordinate2DMake(lat, lon);
}

開發時所面臨的現狀

獲取經緯度(GPS)

  • 火星座標

    • MKMapView

  • 地球座標

    • CLLocationManager

顯示經緯度(地圖)

  • 火星座標

    • iOS 地圖

    • Gogole地圖

    • 搜搜、阿里雲、高德地圖

  • 地球座標

    • Google 衛星地圖(國外地圖應該都是……)

  • 百度座標

    • 百度地圖

推薦的解決方案:

  • 既然是在國內,儲存一律用火星座標,這樣在使用國內地圖顯示時最方便(用百度地圖顯示時可以一次轉換取得)

  • CLLocationManager 拿到的 CLLocation 轉為火星座標,MKMapView 不用處理

  • 使用地圖 API 進行 地址解析/逆地址解析(Geocoding) 時注意相應使用相應地圖商的座標系

  • 部分地圖商支援多個座標系輸入,如高德支援地球、火星座標(這個一直有變動,具體只能參考廠商最新文件了