此转换算法主要用于将国际标准的WGS-84坐标转换为符合中国法律法规的GCJ-02坐标,以便在中国地图服务(如高德、腾讯地图)上正确显示位置。
计算原理分析
1. 椭球体参数定义
a = 6378245.0 # 克拉索夫斯基椭球体长半轴
f = 1 / 298.3 # 椭球体扁率
b = a * (1 - f) # 短半轴
ee = 1 - (b * b)/(a*a) # 第一偏心率平方
这些参数定义了GCJ-02坐标系使用的参考椭球体(克拉索夫斯基椭球体),不同于WGS-84使用的椭球体。
2. 中国境内判断
def outOfChina(lng, lat):
return not (72.004 <= lng <= 137.8347 and 0.8293 <= lat <= 55.8271)
此函数检查坐标是否在中国范围内,只有中国境内的坐标需要转换。
3. 经纬度偏移计算
transformLat
和transformLon
函数使用复杂的多项式和谐波函数计算基础偏移量:
- 包含线性项、二次项和交叉项
- 使用正弦函数模拟周期性偏差
- 算法中的经验公式通过大量实测数据拟合得到
4. 坐标转换核心算法
def wgs2gcj(wgsLon, wgsLat):
if outOfChina(wgsLon, wgsLat):
return wgsLon, wgsLat
# 计算基础偏移量(以105°E, 35°N为参考点)
dLat = transformLat(wgsLon - 105.0, wgsLat - 35.0)
dLon = transformLon(wgsLon - 105.0, wgsLat - 35.0)
# 考虑椭球体形状的修正
radLat = wgsLat / 180.0 * PI
magic = sin(radLat)
magic = 1 - ee * magic * magic
sqrtMagic = sqrt(magic)
# 将偏移量转换为角度
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI)
dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * PI)
# 应用偏移
gcjLat = wgsLat + dLat
gcjLon = wgsLon + dLon
return (gcjLon, gcjLat)
算法特点
- 区域性:只对中国境内的坐标进行转换
- 非线性:偏移量随地理位置变化而变化
- 保密性:使用经验公式而非简单线性变换
- 椭球体修正:考虑了不同参考椭球体之间的差异
需要注意的是,这是一个单向加密过程,精确的逆变换需要更复杂的算法。
WGS84转GCJ-02代码如下:
# -*- coding: utf-8 -*-
from math import sin, cos, sqrt, fabs, atan2
from math import pi as PI
# define ellipsoid
a = 6378245.0
f = 1 / 298.3
b = a * (1 - f)
ee = 1 - (b * b) / (a * a)
def outOfChina(lng, lat):
"""check weather lng and lat out of china
Arguments:
lng {float} -- longitude
lat {float} -- latitude
Returns:
Bollen -- True or False
"""
return not (72.004 <= lng <= 137.8347 and 0.8293 <= lat <= 55.8271)
def transformLat(x, y):
ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(fabs(x))
ret = ret + (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0
ret = ret + (20.0 * sin(y * PI) + 40.0 * sin(y / 3.0 * PI)) * 2.0 / 3.0
ret = ret + (160.0 * sin(y / 12.0 * PI) + 320.0 * sin(y * PI / 30.0)) * 2.0 / 3.0
return ret
def transformLon(x, y):
ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(fabs(x))
ret = ret + (20.0 * sin(6.0 * x * PI) + 20.0 * sin(2.0 * x * PI)) * 2.0 / 3.0
ret = ret + (20.0 * sin(x * PI) + 40.0 * sin(x / 3.0 * PI)) * 2.0 / 3.0
ret = ret + (150.0 * sin(x / 12.0 * PI) + 300.0 * sin(x * PI / 30.0)) * 2.0 / 3.0
return ret
def wgs2gcj(wgsLon, wgsLat):
"""wgs coord to gcj
Arguments:
wgsLon {float} -- lon
wgsLat {float} -- lat
Returns:
tuple -- gcj coords
"""
if outOfChina(wgsLon, wgsLat):
return wgsLon, wgsLat
dLat = transformLat(wgsLon - 105.0, wgsLat - 35.0)
dLon = transformLon(wgsLon - 105.0, wgsLat - 35.0)
radLat = wgsLat / 180.0 * PI
magic = sin(radLat)
magic = 1 - ee * magic * magic
sqrtMagic = sqrt(magic)
dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI)
dLon = (dLon * 180.0) / (a / sqrtMagic * cos(radLat) * PI)
gcjLat = wgsLat + dLat
gcjLon = wgsLon + dLon
return (gcjLon, gcjLat)
点击这里,访问原始代码
更多算法技术问题,欢迎留言或联系我们。转载须注明出处。