图吧地图坐标转GCJ-02（火星）坐标方法



图吧地图坐标转GCJ-02（火星）坐标方法

    

图吧地图坐标转GCJ-02（火星）坐标需要分两步进行：

1）图吧坐标->WGS84坐标。

2）WGS84坐标->火星坐标。这一步在以前是不容易做到的，因为它的算法是一个"秘密"。幸运的是，它已经被某个数学天才"破解"了，在互联网上可以轻易找到其实现算法。

 

具体算法：

1）下面是 图吧坐标转WGS84坐标的JavaScript实现，来源于图吧官方的JS文件。

view plain copy to clipboard print ?

1. function mapBar2WGS84(x, y) {  
2.     x = parseFloat(x) * 100000 % 36000000;  
3.     y = parseFloat(y) * 100000 % 36000000;  
4.   
5.     x1 = parseInt( - (((Math.cos(y / 100000)) * (x / 18000)) + ((Math.sin(x / 100000)) * (y / 9000))) + x);  
6.     y1 = parseInt( - (((Math.sin(y / 100000)) * (x / 18000)) + ((Math.cos(x / 100000)) * (y / 9000))) + y);  
7.   
8.     x2 = parseInt( - (((Math.cos(y1 / 100000)) * (x1 / 18000)) + ((Math.sin(x1 / 100000)) * (y1 / 9000))) + x + ((x > 0) ? 1 : -1));  
9.     y2 = parseInt( - (((Math.sin(y1 / 100000)) * (x1 / 18000)) + ((Math.cos(x1 / 100000)) * (y1 / 9000))) + y + ((y > 0) ? 1 : -1));  
10.   
11.     return [x2 / 100000.0, y2 / 100000.0];  
12. }  

function mapBar2WGS84(x, y) {
    x = parseFloat(x) * 100000 % 36000000;
    y = parseFloat(y) * 100000 % 36000000;

    x1 = parseInt( - (((Math.cos(y / 100000)) * (x / 18000)) + ((Math.sin(x / 100000)) * (y / 9000))) + x);
    y1 = parseInt( - (((Math.sin(y / 100000)) * (x / 18000)) + ((Math.cos(x / 100000)) * (y / 9000))) + y);

    x2 = parseInt( - (((Math.cos(y1 / 100000)) * (x1 / 18000)) + ((Math.sin(x1 / 100000)) * (y1 / 9000))) + x + ((x > 0) ? 1 : -1));
    y2 = parseInt( - (((Math.sin(y1 / 100000)) * (x1 / 18000)) + ((Math.cos(x1 / 100000)) * (y1 / 9000))) + y + ((y > 0) ? 1 : -1));

    return [x2 / 100000.0, y2 / 100000.0];
}

对应的Python实现如下：

view plain copy to clipboard print ?

1. from __future__ import division  
2. import math  
3.    
4. def mapBar2WGS84(lng, lat):  
5.     """图吧坐标转WGS84坐标 
6.     参数： 
7.     lng - 经度 
8.     lat - 纬度 
9.     返回值： 
10.     (经度, 纬度) 
11.     """  
12.     lng = float(lng) * 100000 % 36000000;  
13.     lat = float(lat) * 100000 % 36000000;  
14.     lng1 = int( - (((math.cos(lat / 100000)) * (lng / 18000)) + ((math.sin(lng / 100000)) * (lat / 9000))) + lng)  
15.     lat1 = int( - (((math.sin(lat / 100000)) * (lng / 18000)) + ((math.cos(lng / 100000)) * (lat / 9000))) + lat)  
16.     lng2 = int( - (((math.cos(lat1 / 100000)) * (lng1 / 18000)) + ((math.sin(lng1 / 100000)) * (lat1 / 9000))) + lng + (1 if lng > 0 else -1))  
17.     lat2 = int( - (((math.sin(lat1 / 100000)) * (lng1 / 18000)) + ((math.cos(lng1 / 100000)) * (lat1 / 9000))) + lat + (1 if lat >0 else -1))  
18.     return lng2 / 100000.0, lat2 / 100000.0  

from __future__ import division
import math
 
def mapBar2WGS84(lng, lat):
    """图吧坐标转WGS84坐标
    参数：
    lng - 经度
    lat - 纬度
    返回值：
    (经度, 纬度)
    """
    lng = float(lng) * 100000 % 36000000;
    lat = float(lat) * 100000 % 36000000;
    lng1 = int( - (((math.cos(lat / 100000)) * (lng / 18000)) + ((math.sin(lng / 100000)) * (lat / 9000))) + lng)
    lat1 = int( - (((math.sin(lat / 100000)) * (lng / 18000)) + ((math.cos(lng / 100000)) * (lat / 9000))) + lat)
    lng2 = int( - (((math.cos(lat1 / 100000)) * (lng1 / 18000)) + ((math.sin(lng1 / 100000)) * (lat1 / 9000))) + lng + (1 if lng > 0 else -1))
    lat2 = int( - (((math.sin(lat1 / 100000)) * (lng1 / 18000)) + ((math.cos(lng1 / 100000)) * (lat1 / 9000))) + lat + (1 if lat >0 else -1))
    return lng2 / 100000.0, lat2 / 100000.0

2）WGS84坐标转火星坐标的算法最早成形的版本是 https://on4wp7.codeplex.com/SourceControl/changeset/view/21483#353936，它的Python实现如下：

view plain copy to clipboard print ?

1. from __future__ import division  
2. from math import pi,sqrt,sin,cos  
3.    
4. a = 6378245.0  
5. ee = 0.00669342162296594323  
6.    
7. # World Geodetic System ==> Mars Geodetic System  
8. def transform(wgLat, wgLon):  
9.     """ 
10.     transform(latitude,longitude) , WGS84 
11.     return (latitude,longitude) , GCJ02 
12.     """  
13.     if (outOfChina(wgLat, wgLon)):  
14.         mgLat = wgLat  
15.         mgLon = wgLon  
16.         return  
17.     dLat = transformLat(wgLon - 105.0, wgLat - 35.0)  
18.     dLon = transformLon(wgLon - 105.0, wgLat - 35.0)  
19.     radLat = wgLat / 180.0 * pi  
20.     magic = sin(radLat)  
21.     magic = 1 - ee * magic * magic  
22.     sqrtMagic = sqrt(magic)  
23.     dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi)  
24.     dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi)  
25.     mgLat = wgLat + dLat  
26.     mgLon = wgLon + dLon  
27.     return mgLat,mgLon  
28.    
29. def outOfChina(lat, lon):  
30.     if (lon < 72.004 or lon > 137.8347):  
31.         return True  
32.     if (lat < 0.8293 or lat > 55.8271):  
33.         return True  
34.     return False  
35.    
36. def transformLat(x, y):  
37.     ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x))  
38.     ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0  
39.     ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0  
40.     ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0)) * 2.0 / 3.0  
41.     return ret  
42.    
43. def transformLon(x, y):  
44.     ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x))  
45.     ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0  
46.     ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0  
47.     ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0  
48.     return ret  

from __future__ import division
from math import pi,sqrt,sin,cos
 
a = 6378245.0
ee = 0.00669342162296594323
 
# World Geodetic System ==> Mars Geodetic System
def transform(wgLat, wgLon):
    """
    transform(latitude,longitude) , WGS84
    return (latitude,longitude) , GCJ02
    """
    if (outOfChina(wgLat, wgLon)):
        mgLat = wgLat
        mgLon = wgLon
        return
    dLat = transformLat(wgLon - 105.0, wgLat - 35.0)
    dLon = transformLon(wgLon - 105.0, wgLat - 35.0)
    radLat = wgLat / 180.0 * pi
    magic = sin(radLat)
    magic = 1 - ee * magic * magic
    sqrtMagic = sqrt(magic)
    dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi)
    dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * pi)
    mgLat = wgLat + dLat
    mgLon = wgLon + dLon
    return mgLat,mgLon
 
def outOfChina(lat, lon):
    if (lon < 72.004 or lon > 137.8347):
        return True
    if (lat < 0.8293 or lat > 55.8271):
        return True
    return False
 
def transformLat(x, y):
    ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x))
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0
    ret += (20.0 * sin(y * pi) + 40.0 * sin(y / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * sin(y / 12.0 * pi) + 320 * sin(y * pi / 30.0)) * 2.0 / 3.0
    return ret
 
def transformLon(x, y):
    ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x))
    ret += (20.0 * sin(6.0 * x * pi) + 20.0 * sin(2.0 * x * pi)) * 2.0 / 3.0
    ret += (20.0 * sin(x * pi) + 40.0 * sin(x / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * sin(x / 12.0 * pi) + 300.0 * sin(x / 30.0 * pi)) * 2.0 / 3.0
    return ret

做一个测试：

大众点评的地图系统采用的是图吧地图，就拿我家附近的“小肥羊”（ http://www.dianping.com/shop/550517 ）为例，我们测试一下它在图标地图上的标注和转换为火星坐标后在谷歌地图上的标注点十分一致。

1）它的图吧坐标是：'lat': 34.27071, 'lng': 108.98258。在图吧地图的标注位置如下图所示。

2）计算它的WGS84坐标值：mapBar2WGS84(lng=108.98258, lat=34.27071)的结果为 (108.98525, 34.27116)。

3）计算它的GCj-02坐标（ 火星坐标）：transform(wgLat=34.27116, wgLon=108.98525)的结果为(34.26968285889569, 108.99006153253407)。在谷歌地图(ditu.google.cn)中的标注位置如下图所示。

 

结论：转GCj-02坐标火星坐标后的店与图吧原标注点有很小的偏差（10M内），上述算法精确度很高。