Python GPS数据经纬度坐标纠偏【适用谷歌/高德/百度地图和三大坐标系】

一、三大坐标系

• WGS84坐标系：即地球坐标系，国际上通用的坐标系。
• GCJ02坐标系：即火星坐标系，WGS84坐标系经加密后的坐标系。
• BD09坐标系：即百度坐标系，GCJ02坐标系经加密后的坐标系。

二、纠偏代码
之前做公交路线地图显示，需要根据每路公交车站点的经纬度在地图上标注该站点，所以为了方便起见我将每一路车的经纬度作为一个List进行纠偏处理的而不是单个坐标对，这样处理也比较实用。由于数据是从网络上爬取下来的，List里面的每一个元素都是字符串，为了在处理时方便运算，我把他们进行了处理变成数值类型。代码如下。

import json
import urllib
import math

import numpy
from numpy.core import double

x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方


class Geocoding:
    def __init__(self, api_key):
        self.api_key = api_key

    def geocode(self, address):
        """
        利用高德geocoding服务解析地址获取位置坐标
        :param address:需要解析的地址
        :return:
        """
        geocoding = {'s': 'rsv3',
                     'key': self.api_key,
                     'city': '全国',
                     'address': address}
        geocoding = urllib.urlencode(geocoding)
        ret = urllib.urlopen("%s?%s" % ("http://restapi.amap.com/v3/geocode/geo", geocoding))
        if ret.getcode() == 200:
            res = ret.read()
            json_obj = json.loads(res)
            if json_obj['status'] == '1' and int(json_obj['count']) >= 1:
                geocodes = json_obj['geocodes'][0]
                lng = float(geocodes.get('location').split(',')[0])
                lat = float(geocodes.get('location').split(',')[1])
                return [lng, lat]
            else:
                return None
        else:
            return None

def substring(lng,lat):
#转为数值类型
    tmp_x = []
    tmp_y = []

    for n in lng:
        tmp_x.append(float(n))
    lng = tmp_x;
    for m in lat:
        tmp_y.append(float(m))
    lat = tmp_y;
    return lng,lat



def gcj02_to_bd09(lng, lat):
    """
    火星坐标系(GCJ-02)转百度坐标系(BD-09)
    谷歌、高德——>百度
    :param lng:火星坐标经度
    :param lat:火星坐标纬度
    :return:
    """
    resx = []
    resy = []

    lng,lat=substring(lng,lat)
    for i in range(len(lng)):
        z = math.sqrt(lng[i] * lng[i] + lat[i] * lat[i]) + 0.00002 * math.sin(lat[i] * x_pi)
        theta = math.atan2(lat[i], lng[i]) + 0.000003 * math.cos(lng[i] * x_pi)
        bd_lng = z * math.cos(theta) + 0.0065
        bd_lat = z * math.sin(theta) + 0.006
        resx.append(bd_lng)
        resy.append(bd_lat)
    return resx, resy

def bd09_to_gcj02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    resx = []
    resy = []

    bd_lon,bd_lat=substring(bd_lon,bd_lat)

    for i in range(len(bd_lon)):
        x = bd_lon[i] - 0.0065
        y = bd_lat[i] - 0.006
        z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
        theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
        resx.append(z * math.cos(theta))
        resy.append(z * math.sin(theta))
    return resx, resy


def wgs84_to_gcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """

    resx = []
    resy = []

    lng,lat=substring(lng,lat)
    for i in range(len(lng)):
        dlat = _transformlat(lng[i] - 105.0, lat[i] - 35.0)
        dlng = _transformlng(lng[i] - 105.0, lat[i] - 35.0)
        radlat = lat[i] / 180.0 * pi
        magic = math.sin(radlat)
        magic = 1 - ee * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
        dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
        mglat = lat[i] + dlat
        mglng = lng[i] + dlng
        resx.append(mglat)
        resy.append(mglng)
    return resx, resy




def gcj02_to_wgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    resx = []
    resy = []

    lng,lat=substring(lng,lat)

    for i in range(len(lng)):
        dlat = _transformlat(lng[i] - 105.0, lat[i] - 35.0)
        dlng = _transformlng(lng[i] - 105.0, lat[i] - 35.0)
        radlat = lat[i] / 180.0 * pi
        magic = math.sin(radlat)
        magic = 1 - ee * magic * magic
        sqrtmagic = math.sqrt(magic)
        dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
        dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
        mglat = lat[i] + dlat
        mglng = lng[i] + dlng
        tmpx=lng[i] * 2 - mglng;
        tmpy=lat[i] * 2 - mglat
        resx.append(tmpx)
        resy.append(tmpy)
    return resx, resy


def bd09_to_wgs84(bd_lon, bd_lat):
    lon, lat = bd09_to_gcj02(bd_lon, bd_lat)
    return gcj02_to_wgs84(lon, lat)


def wgs84_to_bd09(lon, lat):
    lon, lat = wgs84_to_gcj02(lon, lat)
    return gcj02_to_bd09(lon, lat)


def _transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
          0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def _transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)


if __name__ == '__main__':
    lngx=['109.485939', '109.48317', '109.480171', '109.478661', '109.478889', '109.479614', '109.481148', '109.483818', '109.484489', '109.484917', '109.488937', '109.490952', '109.492416', '109.493263', '109.491371', '109.491272', '109.487305', '109.482666', '109.478416', '109.476547', '109.472549', '109.466751', '109.461929', '109.461037', '109.461143', '109.456284', '109.450539']
    laty=['36.542042', '36.542622', '36.546299', '36.548508', '36.552353', '36.555672', '36.559906', '36.564976', '36.568703', '36.571568', '36.578037', '36.581493', '36.587299', '36.594879', '36.596882', '36.599083', '36.601543', '36.604279', '36.606976', '36.608147', '36.610722', '36.614956', '36.619328', '36.624397', '36.631588', '36.634911', '36.634281']
    result1 = gcj02_to_bd09(lngx, laty)
    result2 = bd09_to_gcj02(lngx, laty)
    result3 = wgs84_to_gcj02(lngx, laty)
    result4 = gcj02_to_wgs84(lngx, laty)
    result5 = bd09_to_wgs84(lngx, laty)
    result6 = wgs84_to_bd09(lngx, laty)
    lng1,lat1= wgs84_to_gcj02(lngx,laty)
    s,v =gcj02_to_wgs84(lngx,laty)

    print(s)
    print(v)
    print(result1, result2, result3, result4, result5, result6)



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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223