WGS84 与 GCJ-02 坐标系

茶余饭后，我偶然刷到了有关滴滴出行数据的讨论，其中有人提到了 GCJ-02，简单搜索与实践后，将代码与结果整理如下。其中，原始代码参考自 【对错代码分辨】WGS84 坐标系和火星坐标系的转换中的对与错。

WGS84（World Geodetic System 1984）是为 GPS 全球定位系统建立的坐标系统，是世界上第一个统一的地心坐标系，因此也被称为大地坐标系、原始坐标系。一般通过 GPS 记录仪记录下来的经纬度，就是基于 WGS84 坐标系的数据。

GCJ-02 是由中国国家测绘局（G 表示国家，C 表示测绘，J 表示局）制订的地理信息系统的坐标系统，是在 WGS84 经纬度的基础上执行加密算法而成。因为 GPS 得到的经纬度直接在 GCJ-02 坐标系下会定位到错误的地点，有种到了火星的感觉，因此在坊间也将 GCJ-02 戏称为火星坐标系。

WGS84 与 GCJ-02 坐标系间的近似误差（单位：米）
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
import math as m
import typing as t


class Coordinate:
    '''
    - Reference:
        - https://zhuanlan.zhihu.com/p/107253611
    '''

    _a = 6378245.0
    _f = 1 / 298.3
    _b = _a * (1 - _f)
    _ee = 1 - (_b * _b) / (_a * _a)

    def __init__(self, lng: float, lat: float) -> None:
        self._lng, self._lat = lng, lat

    def __repr__(self) -> str:
        return f'<Coordinate-02 @ ({self._lng}, {self._lat})>'

    @classmethod
    def from_gcj(cls, gcj_lng: float, gcj_lat: float) -> 'Coordinate':
        return cls(gcj_lng, gcj_lat)

    @classmethod
    def from_wgs(cls, wgs_lng: float, wgs_lat: float) -> 'Coordinate':
        return cls(*cls._wgs2gcj(wgs_lng, wgs_lat))

    def to_gcj(self) -> t.Tuple[float, float]:
        return self._lng, self._lat

    def to_wgs(self) -> t.Tuple[float, float]:
        return self._gcj2wgs(self._lng, self._lat)

    @classmethod
    def _out_of_china(cls, lng: float, lat: float) -> bool:
        return not (72.004<=lng<=137.8347 and 0.8293<=lat<=55.8271)

    @classmethod
    def _transform_lat(cls, x: float, y: float) -> float:
        return -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*x*y + 0.2*m.sqrt(m.fabs(x)) \
            + (20.0*m.sin(6.0*x*m.pi) + 20.0*m.sin(2.0*x*m.pi))*2.0/3.0 \
            + (20.0*m.sin(y*m.pi) + 40.0*m.sin(y/3.0*m.pi))*2.0/3.0 \
            + (160.0*m.sin(y/12.0*m.pi) + 320.0*m.sin(y*m.pi/30.0))*2.0/3.0

    @classmethod
    def _transform_lng(cls, x: float, y: float) -> float:
        return 300.0 + x + 2.0*y + 0.1*x*x +  0.1*x*y + 0.1*m.sqrt(m.fabs(x)) \
            + (20.0*m.sin(6.0*x*m.pi) + 20.0*m.sin(2.0*x*m.pi))*2.0/3.0 \
            + (20.0*m.sin(x*m.pi) + 40.0*m.sin(x/3.0*m.pi))*2.0/3.0 \
            + (150.0*m.sin(x/12.0*m.pi) + 300.0*m.sin(x*m.pi/30.0))*2.0/3.0

    @classmethod
    def _wgs2gcj(cls, wgs_lng: float, wgs_lat: float) -> t.Tuple[float, float]:
        if cls._out_of_china(wgs_lng, wgs_lat):
            return wgs_lng, wgs_lat
        d_lat = cls._transform_lat(wgs_lng-105.0, wgs_lat-35.0)
        d_lng = cls._transform_lng(wgs_lng-105.0, wgs_lat-35.0)
        rad_lat = wgs_lat / 180.0 * m.pi
        magic = m.sin(rad_lat)
        magic = 1 - cls._ee * magic * magic
        sqrt_magic = m.sqrt(magic)
        d_lat = (d_lat * 180.0) / ((cls._a * (1 - cls._ee)) / (magic * sqrt_magic) * m.pi)
        d_lng = (d_lng * 180.0) / (cls._a / sqrt_magic * m.cos(rad_lat) * m.pi)
        gcj_lat = wgs_lat + d_lat
        gcj_lng = wgs_lng + d_lng
        return gcj_lng, gcj_lat

    @classmethod
    def _gcj2wgs(cls, gcj_lng: float, gcj_lat: float) -> t.Tuple[float, float]:
        g0 = (gcj_lng, gcj_lat)
        w0 = g0
        g1 = cls._wgs2gcj(*w0)
        # w1 = w0 - (g1 - g0)
        w1 = [x[0]-(x[1]-x[2]) for x in zip(w0,g1,g0)]
        # delta = w1 - w0
        delta = [x[0]-x[1] for x in zip(w1, w0)]
        while (abs(delta[0])>=1e-6 or abs(delta[1])>=1e-6):
            w0 = w1
            g1 = cls._wgs2gcj(*w0)
            # w1 = w0 - (g1 - g0)
            w1 = [x[0]-(x[1]-x[2]) for x in zip(w0,g1,g0)]
            # delta = w1 - w0
            delta = [x[0]-x[1] for x in zip(w1, w0)]
        return w1


if __name__ == '__main__':
    wgs_lng, wgs_lat = 112, 40
    gcj_lng, gcj_lat = Coordinate.from_wgs(wgs_lng, wgs_lat).to_gcj()
    wgs_lng2, wgs_lat2 = Coordinate.from_gcj(gcj_lng, gcj_lat).to_wgs()
    delta = m.sqrt((wgs_lng2-wgs_lng)**2 + (wgs_lat2-wgs_lat)**2)
    print(f'{delta = }')

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import itertools as it

import matplotlib.pyplot as plt
import numpy as np

from mpl_toolkits.basemap import Basemap


class Coordinate:
    ...


n_lng, n_lat = 256, 256
wgs_lngs, wgs_lats = np.meshgrid(
    np.linspace(72.004, 137.8347, n_lng), np.linspace(0.8293, 55.8271, n_lat)
)
gcj_lngs, gcj_lats = np.zeros_like(wgs_lngs), np.zeros_like(wgs_lats)
for ith, jth in it.product(range(n_lat), range(n_lng)):
    gcj_lngs[ith, jth], gcj_lats[ith, jth] = \
        Coordinate.from_wgs(wgs_lngs[ith, jth], wgs_lats[ith, jth]).to_gcj()
delta_lng = (wgs_lngs-gcj_lngs) / 1.141255544679108e-5
delta_lat = (wgs_lats-gcj_lats) / 8.993216192195822e-6
delta = np.sqrt(delta_lng**2 + delta_lat**2)

fig, ax = plt.subplots(1, figsize=(16, 9))
# setup mercator map projection.
m = Basemap(
    llcrnrlon=72.004, llcrnrlat=0.8293, urcrnrlon=137.8347, urcrnrlat=55.8271,
    rsphere=(6378137.00, 6356752.3142),
    resolution='l', projection='merc',
    lat_0=40., lon_0=-20., lat_ts=20.,
)
m.drawcoastlines()
m.drawcountries()
# draw parallels and meridians
m.drawparallels(np.arange(10, 90, 10),labels=[1, 1, 0, 1])
m.drawmeridians(np.arange(-180, 180, 10),labels=[1, 1, 0, 1])
im = m.imshow(delta[::-1, :], interpolation='nearest')
fig.colorbar(im, ax=ax)
plt.show()