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