python代码库
库地址：https://pypi.org/project/coord-convert/
源作者Github首页：https://github.com/sshuair/coord-convert
什么开源协助咱也不懂，咱也不是代码原创，源作者地址就挂这里。[b]本贴主要分享思路
首先是安装库，不多介绍
找到核心转换源码
简化python核心代码
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):
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):
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)
def gcj2wgs(gcjLon, gcjLat):
g0 = (gcjLon, gcjLat)
w0 = g0
g1 = wgs2gcj(w0[0], w0[1])
# w1 = w0 - (g1 - g0)
w1 = tuple(map(lambda x: x[0]-(x[1]-x[2]), zip(w0,g1,g0)))
# delta = w1 - w0
delta = tuple(map(lambda x: x[0] - x[1], zip(w1, w0)))
while (abs(delta[0]) >= 1e-6 or abs(delta[1]) >= 1e-6):
w0 = w1
g1 = wgs2gcj(w0[0], w0[1])
# w1 = w0 - (g1 - g0)
w1 = tuple(map(lambda x: x[0]-(x[1]-x[2]), zip(w0,g1,g0)))
# delta = w1 - w0
delta = tuple(map(lambda x: x[0] - x[1], zip(w1, w0)))
return w1
Excel中导入使用
插入py，并粘贴代码进去，即可。
调用上述python代码中的函数
得到的是”元组“
元组展开
展开后转置下
将python核心代码交给ChatGPT，转换成VBA代码
因为from导入的模块，都可以用VBA直接实现，因此不存在什么技术壁垒，ChatGPT可以直接进行转换
from math import sin, cos, sqrt, fabs, atan2
from math import pi as PI
得到VBA代码
作者删除了VBA中部分As Double类型声明，避免因类型不同而报错。
' Constants
Const PI As Double = 3.14159265358979
' Define ellipsoid
Const a As Double = 6378245#
Const f As Double = 1 / 298.3
Const b As Double = a * (1 - f)
Const ee As Double = 1 - (b * b) / (a * a)
Function outOfChina(lng, lat) As Boolean
outOfChina = Not (72.004 <= lng And lng <= 137.8347 And 0.8293 <= lat And lat <= 55.8271)
End Function
Function transformLat(x, y)
Dim ret As Double
ret = -100# + 2# * x + 3# * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Sqr(Abs(x))
ret = ret + (20# * Sin(6# * x * PI) + 20# * Sin(2# * x * PI)) * 2# / 3#
ret = ret + (20# * Sin(y * PI) + 40# * Sin(y / 3# * PI)) * 2# / 3#
ret = ret + (160# * Sin(y / 12# * PI) + 320# * Sin(y * PI / 30#)) * 2# / 3#
transformLat = ret
End Function
Function transformLon(x, y)
Dim ret As Double
ret = 300# + x + 2# * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Sqr(Abs(x))
ret = ret + (20# * Sin(6# * x * PI) + 20# * Sin(2# * x * PI)) * 2# / 3#
ret = ret + (20# * Sin(x * PI) + 40# * Sin(x / 3# * PI)) * 2# / 3#
ret = ret + (150# * Sin(x / 12# * PI) + 300# * Sin(x * PI / 30#)) * 2# / 3#
transformLon = ret
End Function
Function wgs2gcj(wgsLon, wgsLat) As Variant
Dim dLat As Double
Dim dLon As Double
Dim radLat As Double
Dim magic As Double
Dim sqrtMagic As Double
Dim gcjLon As Double
Dim gcjLat As Double
If outOfChina(wgsLon, wgsLat) Then
gcjLon = wgsLon
gcjLat = wgsLat
Else
dLat = transformLat(wgsLon - 105#, wgsLat - 35#)
dLon = transformLon(wgsLon - 105#, wgsLat - 35#)
radLat = wgsLat / 180# * PI
magic = Sin(radLat)
magic = 1 - ee * magic * magic
sqrtMagic = Sqr(magic)
dLat = (dLat * 180#) / ((a * (1 - ee)) / (magic * sqrtMagic) * PI)
dLon = (dLon * 180#) / (a / sqrtMagic * Cos(radLat) * PI)
gcjLat = wgsLat + dLat
gcjLon = wgsLon + dLon
End If
wgs2gcj = Array(gcjLon, gcjLat)
End Function
Function gcj2wgs(gcjLon, gcjLat) As Variant
Dim g0 As Variant
Dim w0 As Variant
Dim g1 As Variant
Dim w1 As Variant
Dim delta As Variant
g0 = Array(gcjLon, gcjLat)
w0 = g0
g1 = wgs2gcj(w0(0), w0(1))
w1 = Array(w0(0) - (g1(0) - g0(0)), w0(1) - (g1(1) - g0(1)))
delta = Array(w1(0) - w0(0), w1(1) - w0(1))
Do While Abs(delta(0)) >= 0.000001 Or Abs(delta(1)) >= 0.000001
w0 = w1
g1 = wgs2gcj(w0(0), w0(1))
w1 = Array(w0(0) - (g1(0) - g0(0)), w0(1) - (g1(1) - g0(1)))
delta = Array(w1(0) - w0(0), w1(1) - w0(1))
Loop
gcj2wgs = w1
End Function
奥维测试
随便找个点GCJ坐标：113.95566478，30.42707824
这个是Excel运行python代码和VBA代码，分别得到的结果，基本一致。
奥维测试结果两者偏差，肉眼不可见。基本认为此代码可用。
结论
我就是想使用VBA代码，直接使用VBA代码就行了，不知道有什么更好的方法，欢迎讨论。
免责声明：测试不严谨，数量也不多，不能保证完全不出错。
发表于 2024-4-15 16:02
