1 为什么会有高德坐标系和百度坐标系？

• 根据《测绘法》和国家保密法规，在中国大陆范围内的地理坐标数据必须做加密处理，不允许直接使用 WGS84（openstreetmap）
  • 所以出现了GCJ-02 和 BD-09
    • 高德、腾讯、谷歌中国都遵循 GCJ-02
      • 在 WGS-84 的基础上，加了随机扰动和非线性偏移算法（保密算法）
      • 火星坐标系，是一种基于WGS-84制定的大地测量系统，由中国国测局制定。
      • 此坐标系所采用的混淆算法会在经纬度中加入随机的偏移。
    • 百度在 GCJ-02 的基础上又加了一次二次加密，形成 BD-09
      • 在 GCJ-02 的基础上，再做一次偏移和加密，主要是 保护商业利益和数据差异化
    • 因此，不同平台互换数据时，需要做坐标系转换，否则会出现“偏移”。
• WGS-84（World Geodetic System, WGS）是使用最广泛的坐标系，也是世界通用的坐标系，GPS设备得到的经纬度就是在WGS84坐标系下的经纬度。

2  百度坐标转高德坐标

import mathdef bd09_to_gcj02(bd_lon, bd_lat):"""百度坐标系 (BD09) 转 高德/谷歌 (GCJ02)"""x = bd_lon - 0.0065y = bd_lat - 0.006z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * math.pi)theta = math.atan2(y, x) - 0.000003 * math.cos(x * math.pi)gg_lon = z * math.cos(theta)gg_lat = z * math.sin(theta)return gg_lon, gg_lat# 示例：百度坐标 -> 高德坐标
bd_lon, bd_lat = 121.495, 31.240   # 上海某点的百度坐标
gcj_lon, gcj_lat = bd09_to_gcj02(bd_lon, bd_lat)print("百度坐标：", bd_lon, bd_lat)
#百度坐标： 121.495 31.24
print("高德坐标：", gcj_lon, gcj_lat)
#高德坐标： 121.48850960667689 31.23401650463438

3  高德地图转百度地图

import mathdef gcj02_to_bd09(gcj_lon, gcj_lat):"""高德/腾讯/谷歌中国 (GCJ-02) → 百度 (BD-09)"""z = math.sqrt(gcj_lon * gcj_lon + gcj_lat * gcj_lat) + 0.00002 * math.sin(gcj_lat * math.pi)theta = math.atan2(gcj_lat, gcj_lon) + 0.000003 * math.cos(gcj_lon * math.pi)bd_lon = z * math.cos(theta) + 0.0065bd_lat = z * math.sin(theta) + 0.006return bd_lon, bd_lat# 示例
gcj_lon, gcj_lat = 121.48, 31.23
bd_lon, bd_lat = gcj02_to_bd09(gcj_lon, gcj_lat)
print(bd_lon, bd_lat)
#121.48649307312748 31.23597382351204

4 高德地图转WGS-84

import mathPI = math.pi
A = 6378245.0  # 长半轴
EE = 0.00669342162296594323  # 偏心率平方def transform_lat(x, y):ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + \0.1 * x * y + 0.2 * math.sqrt(abs(x))ret += (20.0 * math.sin(6.0 * x * PI) + 20.0 *math.sin(2.0 * x * PI)) * 2.0 / 3.0ret += (20.0 * math.sin(y * PI) + 40.0 *math.sin(y / 3.0 * PI)) * 2.0 / 3.0ret += (160.0 * math.sin(y / 12.0 * PI) + 320 *math.sin(y * PI / 30.0)) * 2.0 / 3.0return retdef transform_lon(x, y):ret = 300.0 + x + 2.0 * y + 0.1 * x * x + \0.1 * x * y + 0.1 * math.sqrt(abs(x))ret += (20.0 * math.sin(6.0 * x * PI) + 20.0 *math.sin(2.0 * x * PI)) * 2.0 / 3.0ret += (20.0 * math.sin(x * PI) + 40.0 *math.sin(x / 3.0 * PI)) * 2.0 / 3.0ret += (150.0 * math.sin(x / 12.0 * PI) + 300.0 *math.sin(x / 30.0 * PI)) * 2.0 / 3.0return retdef gcj02_to_wgs84(lon, lat):dlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * PImagic = math.sin(radlat)magic = 1 - EE * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * PI)dlon = (dlon * 180.0) / (A / sqrtmagic * math.cos(radlat) * PI)mgLat = lat + dlatmgLon = lon + dlonreturn lon * 2 - mgLon, lat * 2 - mgLatgcj_lon, gcj_lat = 121.48, 31.23  # 高德
wgs_lon, wgs_lat = gcj02_to_wgs84(gcj_lon, gcj_lat)
print("高德(GCJ02):", gcj_lon, gcj_lat)
print("OSM(WGS84):", wgs_lon, wgs_lat)
'''
高德(GCJ02): 121.48 31.23
OSM(WGS84): 121.47549797745022 31.231959998113442
'''

5 WGS84 转高德

import mathdef transform_lat(lon, lat):ret = -100.0 + 2.0 * lon + 3.0 * lat + 0.2 * lat * lat + \0.1 * lon * lat + 0.2 * math.sqrt(abs(lon))ret += (20.0 * math.sin(6.0 * lon * math.pi) + 20.0 * math.sin(2.0 * lon * math.pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lat * math.pi) + 40.0 * math.sin(lat / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(lat / 12.0 * math.pi) + 320.0 * math.sin(lat * math.pi / 30.0)) * 2.0 / 3.0return retdef transform_lon(lon, lat):ret = 300.0 + lon + 2.0 * lat + 0.1 * lon * lon + \0.1 * lon * lat + 0.1 * math.sqrt(abs(lon))ret += (20.0 * math.sin(6.0 * lon * math.pi) + 20.0 * math.sin(2.0 * lon * math.pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lon * math.pi) + 40.0 * math.sin(lon / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(lon / 12.0 * math.pi) + 300.0 * math.sin(lon / 30.0 * math.pi)) * 2.0 / 3.0return retdef wgs84_to_gcj02(lat,lon):"""WGS84 → GCJ02"""dlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - 0.00669342162296594323 * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((6378245.0 * (1 - 0.00669342162296594323)) / (magic * sqrtmagic) * math.pi)dlon = (dlon * 180.0) / (6378245.0 / sqrtmagic * math.cos(radlat) * math.pi)mglat = lat + dlatmglon = lon + dlonreturn mglat,mglon