GCJ坐标转换及其生成VBA代码思路

小能维尼

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代码就行了，不知道有什么更好的方法，欢迎讨论。

免责声明：测试不严谨，数量也不多，不能保证完全不出错。

niccccce

受教了！！