原始GPS坐标 (WGS-84)转换为 GCJ02坐标 和 BD09坐标 python源码

本文地址：http://tongxinmao.com/Article/Detail/id/510

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
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
# coding: utf-8
from __future__ import division
import re
import json
import urllib
import math
 
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 gcj02_to_bd09(lng, lat):
    """
    火星坐标系(GCJ-02)转百度坐标系(BD-09)
    谷歌、高德——>百度
    :param lng:火星坐标经度
    :param lat:火星坐标纬度
    :return:
    """
    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)
    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)
    bd_lng = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    return [bd_lng, bd_lat]
 
 
def bd09_to_gcj02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    x = bd_lon - 0.0065
    y = bd_lat - 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)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return [gg_lng, gg_lat]
 
 
def wgs84_to_gcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
 
 
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return lng, lat
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 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 + dlat
    mglng = lng + dlng
    return [mglng, mglat]
 
 
def gcj02_to_wgs84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return lng, lat
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 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 + dlat
    mglng = lng + dlng
    return [lng * 2 - mglng, lat * 2 - mglat]
 
 
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)
     
def ddmmsstodd(lng, lat):
    #度分秒转度（带°′"符号，如120°25′17"）
    lng = re.split("°|′|″", lng)[:3]
    lat = re.split("°|′|″", lat)[:3]
    olng = list(map(int, lng))
    olat = list(map(int, lat))
    return [(olng[0] + olng[1] / 60 + olng[2] / 3600), (olat[0] + olat[1] / 60 + olat[2] / 3600)]
     
def ddmmtoddd(lng, lat):
    lng = str(lng/100)
    lng = lng.split(".")[0] + "." + str(int((lng.split("."))[1])/100*100/60).replace(".", "")[:9]
     
    lat = str(lat/100)
    lat = lat.split(".")[0] + "." + str(int((lat.split("."))[1])/100*100/60).replace(".", "")[:9]
    return [lng, lat]
 
 
if __name__ == '__main__':
 
    print("")
    print("coordinate convert:")
    print("-----------------")
 
    #dd.dddd
    lng = 128.543
    lat = 22.5635350
    print("lng", lng)
    print("lat", lat)
    print("GCJ-02 to BD-09", gcj02_to_bd09(lng, lat))       #火星坐标系(GCJ-02)转百度坐标系(BD-09)
    print("BD-09 to GCJ-02", bd09_to_gcj02(lng, lat))       #百度坐标系(BD-09)转火星坐标系(GCJ-02)
    print("WGS-84 to GCJ-02", wgs84_to_gcj02(lng, lat))     #WGS84转GCJ02(火星坐标系)
    print("GCJ-02 to WGS-84", gcj02_to_wgs84(lng, lat))     #GCJ02(火星坐标系)转GPS84
    print("BD-09 to WGS-84", bd09_to_wgs84(lng, lat))       #百度坐标系(BD-09)转GPS84
    print("WGS-84 to BD-09", wgs84_to_bd09(lng, lat))       #WGS84转百度坐标系(BD-09)
 
    #g = Geocoding('API_KEY')  # 这里填写你的高德api的key
     
    #dd°mm′ss.ss″ to dd.dddd
    olng = "120°25′17″"
    olat = "22°18′81″"
    print("")
    print("dd°mm′ss.ss″ to dd.dddddd")
    print("-----------------")
    print("olng", olng)
    print("olat", olat)
    print(ddmmsstodd(olng, olat))
     
     
    #dd.dddd to dd°mm′ss.ss″
    #to do
     
    #ddmm.mmmmm to dd.dddddd
    lng = 11355.8748
    lat = 2233.8121
     
    '''unction dm2dd(f: double): double;
        begin
          Result := Int(f / 100) + Frac(f / 100) * 100 / 60;
        end;'''
    print("")
    print("ddmm.mmmm″ to dd.dddddd")
    print("-----------------") 
    print(lng, lat)
    print(ddmmtoddd(lng, lat))
     
     
    print("")
    print("ddmm.mmmm″ to dd.dddddd to BD-09")
    print("-----------------") 
    lng = 11355.8748
    lat = 2233.8121
    print()
    print("转百度坐标系：", wgs84_to_bd09(float(ddmmtoddd(lng, lat)[0]), float(ddmmtoddd(lng, lat)[1])))
    print("")
    print("")
    print("")
     
    inputs= input("请输入lng, lat （ddmm.mmmm格式）, 以 ',' 分隔：")
    if "," in inputs:
        print("ddmm.mmmm 转 dd.dddd")
        print(ddmmtoddd(float(inputs.split(",")[0]), float(inputs.split(",")[1])))
        alat = float(ddmmtoddd(float(inputs.split(",")[0]), float(inputs.split(",")[1]))[0])
        alng = float(ddmmtoddd(float(inputs.split(",")[0]), float(inputs.split(",")[1]))[1])
        print(alng,alat)
         
        print("")
        print("转百度坐标系：", wgs84_to_bd09(alng, alat))
    else:
        print("输入错误！")


上一篇：利用虚拟桌面无窗口启动进程
下一篇：LUNA: a USB multitool & nMigen library