坐标转换与地理编码

坐标系统,是描述物质存在的空间位置(坐标)的参照系,通过定义特定基准及其参数形式来实现。坐标是描述位置的一组数值,按坐标的维度一般分为一维坐标(公路里程碑)和二维坐标(笛卡尔平面直角坐标)、三维坐标(大地坐标、空间直角坐标)。为了描述或确定位置,必须建立坐标系统,坐标只有存在于某个坐标系统才有实际的意义与具体的位置。
地理坐标系统:地理坐标坐标是对地球进行简单几何建模,比如将地球看成一个球体或者类球体,然后再将地表上点投影到该球面上形成的坐标就是地理坐标系统。
投影坐标系统:由于地球是一个球状,所以一般将其某个区域投影在平面上,形成的坐标系称为投影坐标系。

坐标转换

WGS-84坐标系:地理坐标系统,即GPS原始坐标体系。目前基本上所有定位空间位置的设备都使用这种坐标系统,例如手机的GPS系统。
GCJ-02坐标系:投影坐标系统,也就是我们平常所说的火星坐标系,这个是中国自己在WGS84基础上加密而成。
BD-09坐标系:投影坐标系统,百度在GCJ-02的基础上,又做了一次加密,又称BD-09坐标系。

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
import math

#========================================================
# 坐标转换
#========================================================
a = 6378245.0 # 克拉索夫斯基椭球参数长半轴a
f = 1/293.3
ee = 2 * f - f * f # 克拉索夫斯基椭球参数第一偏心率平方
bd_pi = math.pi * 3000.0 / 180.0

def standard_loc(loc):
'''
度分秒坐标转化为度
INPUT -> 字符串格式坐标(度分秒)
OUTPUT -> 坐标(度)
'''
loc = str(loc)
degree = float(loc[:loc.find('°')])
minute = float(loc[loc.find('°') + 1:loc.find('′')])/60
second = float(loc[loc.find('′') + 1:loc.find('″')])/3600
return degree+minute+second

def out_of_china(wgsLon, wgsLat):
'''
判断是否在国内
INPUT -> GPS经度, GPS纬度
OUTPUT -> 是-不在国内/否-在国内
'''
return not (wgsLon > 72.004 and wgsLon < 137.8347 and wgsLat > 0.8293 and wgsLat < 55.8271 )

def gcj_transformLon(x, y):
'''
转换GCJ经度偏移量
'''
ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * y * x + 0.1 * math.sqrt(math.fabs(x))
ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
return ret

def gcj_transformLat(x, y):
'''
计算GCJ纬度偏移量
'''
ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * y * x + 0.2 * math.sqrt(math.fabs(x))
ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
return ret

def gcj_delta(lon, lat):
'''
计算GCJ坐标偏移量
'''
radLat = lat * math.pi / 180.0
magic = 1 - ee * math.sin(radLat) * math.sin(radLat)
sqrtMagic = math.sqrt(magic)
dLon = gcj_transformLon(lon - 105.0, lat - 35.0)
dLat = gcj_transformLat(lon - 105.0, lat - 35.0)
dLon = (dLon * 180.0) / (a / sqrtMagic * math.cos(radLat) * math.pi)
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * math.pi)
return dLon, dLat

def wgs84_to_gcj02(wgsLon, wgsLat):
'''
WGS84(地球坐标)转GCJ02(火星坐标系)
INPUT -> GPS经度, GPS纬度
OUTPUT -> GCJ经度, GCJ纬度
'''
if out_of_china(wgsLon, wgsLat): # 判断是否在国外
print(u'该地点不在中国境内,不执行转换!')
return wgsLon, wgsLat
dLon, dLat = gcj_delta(wgsLon, wgsLat)
# 得到偏移坐标
gcjLon = wgsLon + dLon
gcjLat = wgsLat + dLat
return gcjLon, gcjLat

def gcj02_to_bd09(gcjLon, gcjLat):
'''
GCJ02(火星坐标系)转BD09(百度坐标)
INPUT -> GCJ经度, GCJ纬度
OUTPUT -> BD经度, BD纬度
'''
x = gcjLon
y = gcjLat
z = math.sqrt(x * x + y * y) + 0.00002 * math.sin(y * bd_pi)
theta = math.atan2(y, x) + 0.000003 * math.cos(x * bd_pi)
bdLon = z * math.cos(theta) + 0.0065
bdLat = z * math.sin(theta) + 0.006
return bdLon, bdLat

def bd09_to_gcj02(bdlon, bdLat):
'''
BD09(百度坐标转)GCJ02(火星坐标系)
INPUT -> BD经度, BD纬度
OUTPUT -> GCJ经度, GCJ纬度
'''
x = bdlon - 0.0065
y = bdLat - 0.006
z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * bd_pi)
theta = math.atan2(y, x) - 0.000003 * math.cos(x * bd_pi)
gcjLon = z * math.cos(theta)
gcjLat = z * math.sin(theta)
return gcjLon, gcjLat

def gcj02_to_wgs84(gcjLon, gcjLat):
'''
GCJ02(火星坐标系)转WGS84(地球坐标)
INPUT -> GCJ经度, GCJ纬度
OUTPUT -> GPS经度, GPS纬度
'''
dLon, dLat = gcj_delta(gcjLon, gcjLat)
# 得到偏移坐标
gcjLon = gcjLon - dLon
gcjLat = gcjLat - dLat
return gcjLon, gcjLat

def dis_calculate(x1, y1, x2, y2):
'''
球面距离计算
INPUT -> 坐标1经度, 坐标1纬度, 坐标2经度, 坐标2纬度
OUTPUT -> 两点间球面距离
'''
return math.sqrt(math.pow(x1-x2,2)+math.pow(y1-y2,2))*3600*26.5

地理编码

我们应用电子地图查询信息,有时候希望获取点位坐标,有时候希望得到一个准确的地址,这在地图底层服务里面实际是用到了地理编码和逆地理编码。
地理编码是指将地址或地名等位置描述转换为经纬度坐标的过程。得到的坐标信息,可以用于制图或空间分析操作。通过地理编码可快速查找到各类位置。比如建筑物名称、地址、位置描述、名胜景区等等。
逆地理编码实际是一个空间查询的过程,通过输入经纬度坐标,查找这个坐标所在的行政区划,所属道路,最近的POI/AOI等。这些数据一般比较庞杂,所以合理组织这些数据结构非常关键。管理空间位置数据一般用RTree结构。RTree是一种高度平衡的树,在叶子节点存放数据指针。使用RTree能保证对一个空间数据的搜素只需要访问很小一部分节点就能查找到数据。

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
import requests
import json
#========================================================
# 地理编码与逆地理编码
#========================================================

def getlonlat_by_address(address):
'''
地理编码:将地理名词通过高德地图转换为经纬度
INPUT -> 地理名词
OUTPUT -> 高德坐标
'''
url = 'https://restapi.amap.com/v3/geocode/geo'
key = '5bdc284755f32a8f48a5c179507a551d'
uri = url+'?'+'address='+address+'&output=json&key='+key
res = requests.get(uri).text
result = json.loads(res)['geocodes'][0]['location']
return result

def getaddress_by_lonlat(lonlat):
'''
逆地理编码:将经纬度通过高德地图转换为地理名词
INPUT -> 高德坐标
OUTPUT -> 地理名词
'''
url = 'https://restapi.amap.com/v3/geocode/regeo'
key = '5bdc284755f32a8f48a5c179507a551d'
uri = url+'?'+'location='+lonlat+'&output=json&radius=1000&extensions=all&key='+key
res = requests.get(uri).text
result = json.loads(res)
return result
0%