[VBA代码]高斯正反算（经纬度和平面坐标转换）

小能维尼

[VBA]高斯正反算

代码来源

利用Python编写一个高斯正反算程序-CSDN博客

声明：本人仅是"简单"修改和润色（和chatgpt帮助）

分享给各位需要的朋友

python的库代码看不懂，又找不到核心代码
网上找到其他的python代码计算结果不对，难搞

VBA代码

Dim a As Double, e1 As Double, e2 As Double
Dim N1 As Integer, L0 As Double, L1 As Double

Const pi As Double = 3.14159265358979

' 固定参数
Const 椭球参数str As String = "CGCS2000"
Const 带号选择str As String = "三度带"

' 选择椭球参数
Private Function EllipsoidParameter(str1 As String) As Variant
    Dim m0 As Double, m2 As Double, m4 As Double, m6 As Double, m8 As Double
    Dim a0 As Double, a2 As Double, a4 As Double, a6 As Double, a8 As Double

    If str1 = "CGCS2000" Then
        a = 6378137
        B = 6356752.31414036
    ElseIf str1 = "WGS84" Then
        a = 6378137
        B = 6356752.31424518
    ElseIf str1 = "1980国家大地坐标系" Then
        a = 6378140
        B = 6356755.28815753
    ElseIf str1 = "克拉索夫斯基椭球体" Then
        a = 6378245
        B = 6356863.01877305
    End If

    e1 = (a * a - B * B) / (a * a)
    e2 = (a * a - B * B) / (B * B)

    m0 = a * (1 - e1)
    m2 = 1.5 * m0 * e1
    m4 = 1.25 * m2 * e1
    m6 = (7# / 6) * m4 * e1
    m8 = (9# / 8) * m6 * e1

    a0 = m0 + m2 / 2 + (3# / 8) * m4 + (5# / 16) * m6 + (35# / 128) * m8
    a2 = (1# / 2) * (m2 + m4) + (15# / 32) * m6 + (7# / 16) * m8
    a4 = (1# / 8) * m4 + (3# / 16) * m6 + (7# / 32) * m8
    a6 = (1# / 32) * m6 + (1# / 16) * m8
    a8 = (1# / 128) * m8

    EllipsoidParameter = Array(a, e1, e2, a0, a2, a4, a6, a8)
End Function

' 大地转高斯带号选择
Private Function NumberedSelectionBL2XY(str2 As String, L As Double) As Variant
    If str2 = "三度带" Then
        N1 = Int((L + 1.5) / 3)
        L0 = 3 * N1
    ElseIf str2 = "六度带" Then
        N1 = Int(L / 6) + 1
        L0 = 6 * N1 - 3
    End If
    NumberedSelectionBL2XY = Array(N1, L0)  '带号 和 中央经线
End Function

' 高斯转大地带号选择
Private Function NumberedSelectionXY2BL(str2 As String, aa As Long) As Double
    If str2 = "三度带" Then
        L1 = 3 * aa
    ElseIf str2 = "六度带" Then
        L1 = 6 * aa - 3
    End If
    NumberedSelectionXY2BL = L1
End Function

' 输入十进制度 经度 纬度
Function LB2XY(L As Double, B As Double, Optional DH As Integer = -1) As Variant
    Dim list1 As Variant
    Dim list2 As Variant
    Dim radB As Double, delta1 As Double, delta2 As Double
    Dim delta3 As Double, delta4 As Double, delta5 As Double
    Dim x As Double, n As Double, t As Double, yita As Double
    Dim roum As Double, l_1 As Double, x_1 As Double, y As Double

    ' 选择椭球参数
    list1 = EllipsoidParameter(椭球参数str)

    ' 大地转高斯带号选择
    list2 = NumberedSelectionBL2XY(带号选择str, L)

    ' 强制带号转换
    If DH <> -1 Then
        list2(0) = DH           '强制-带号
        list2(1) = 3 * DH       '强制-中央经线
    End If

    ' 求椭球面上某点到赤道的子午线弧长，精度至0.001m
    radB = (pi / 180) * B
    delta1 = list1(3) * radB
    delta2 = (list1(4) / 2) * Sin(2 * radB)
    delta3 = (list1(5) / 4) * Sin(4 * radB)
    delta4 = (list1(6) / 6) * Sin(6 * radB)
    delta5 = (list1(7) / 8) * Sin(8 * radB)
    x_1 = Round(delta1 - delta2 + delta3 - delta4 + delta5, 9)

    n = list1(0) / Sqr(1 - list1(1) * Sin(radB) * Sin(radB))
    t = Tan(radB)
    yita = Sqr(list1(2)) * Cos(radB)
    roum = 180# * 60# * 60# / pi
    l_1 = ((L - list2(1)) * 3600) / roum

    x = Round(x_1 + (n / 2) * t * Cos(radB) * Cos(radB) * l_1 * l_1 + _
              (n / 24) * t * (5 - t * t + 9 * yita * yita + 4 * yita ^ 4) * Cos(radB) ^ 4 * l_1 ^ 4 + _
              (n / 720) * t * (61 - 58 * t * t + t ^ 4) * l_1 ^ 6 * Cos(radB) ^ 6, 4)

    y = Round(n * Cos(radB) * l_1 + _
              (n / 6) * (1 - t * t + yita * yita) * Cos(radB) ^ 3 * l_1 ^ 3 + _
              (n / 120) * (5 - 18 * t * t + t ^ 4 + 14 * yita * yita - 58 * yita * yita * t * t) * _
              Cos(radB) ^ 5 * l_1 ^ 5 + 500000 + list2(0) * 10 ^ 6, 4)
    '地理坐标系
    LB2XY = Array(x, y)
End Function

' 高斯坐标转大地坐标
Function XY2LB(x As Double, y1 As Double) As Variant

    Dim list1 As Variant
    Dim y As Double, aa As Long
    Dim Bf1(19) As Double, F(19) As Double
    Dim deltad As Double, tf As Double, yitaf As Double, Mf As Double, Nf As Double
    Dim i As Integer
    Dim L As Double, degrees1 As Double, B As Double, degrees2 As Double
    Dim L1 As Double

    ' 获取CGCS2000的椭球参数
    list1 = EllipsoidParameter(椭球参数str)

    ' 计算y的值
    y = y1 - Int(y1 / 10 ^ 6) * 10 ^ 6 - 500000
    aa = Int(y1 / 10 ^ 6)

    ' 底点纬度
    Bf1(0) = x / list1(3)

    ' 迭代求解Bf1
    For i = 0 To 9
        F(i) = (-list1(4) / 2) * Sin(2 * Bf1(i)) + list1(5) / 4 * Sin(4 * Bf1(i)) - list1(6) / 6 * Sin(6 * Bf1(i))
        Bf1(i + 1) = (x - F(i)) / list1(3)
        deltad = Bf1(i + 1) - Bf1(i)
        If Abs(deltad) < 10 ^ -10 Then
            Exit For
        End If
    Next i

    ' 计算高斯坐标系参数
    tf = Tan(Bf1(i))
    yitaf = list1(2) * Cos(Bf1(i)) ^ 2
    Mf = list1(0) * (1 - list1(1)) / (1 - list1(1) * Sin(Bf1(i)) ^ 2) ^ 1.5
    Nf = list1(0) / Sqr(1 - list1(1) * Sin(Bf1(i)) ^ 2)

    ' 带号选择
    L1 = NumberedSelectionXY2BL(带号选择str, aa)

    ' 计算经度L和纬度B
    L = (WorksheetFunction.pi() / 180) * L1 + y / (Nf * Cos(Bf1(i))) - (1 + 2 * tf ^ 2 + yitaf) * y ^ 3 / (6 * Nf ^ 3 * Cos(Bf1(i))) + (5 + 28 * tf ^ 2 + 24 * tf ^ 4 + 6 * yitaf + 8 * yitaf * tf ^ 4) * y ^ 5 / (120 * Nf ^ 5 * Cos(Bf1(i)))
    degrees1 = (180 / WorksheetFunction.pi()) * L

    B = Bf1(i) - (tf * y ^ 2) / (2 * Mf * Nf) + tf * (5 + 3 * tf ^ 2 + yitaf - 9 * yitaf * tf ^ 2) * y ^ 4 / (24 * Mf * Nf ^ 3) + tf * (61 + 90 * tf ^ 2 + 45 * tf ^ 4) * y ^ 6 / (720 * Mf * Nf ^ 5)
    degrees2 = (180 / WorksheetFunction.pi()) * B

    ' 返回经纬度
    XY2LB = Array(degrees1, degrees2)
End Function

使用方法

固定参数：在vba代码中，已经固定参数，够用了，懂的可自行修改。

' 固定参数
Const 椭球参数str As String = "CGCS2000"
Const 带号选择str As String = "三度带"

函数：

L，B = XY2LB(x值,y值)

X，Y = LB2XY(经度,纬度,[带号])   

*注0：函数返回数组，依靠Excel动态数组溢出功能显示，旧版本可能需要特殊方法

*注1：带号为可选项，默认根据经度计算带号，可强制设置带号

*注2：经纬度为十进制度的单位

注意：XY结果的计算精度在0.001m

isunfly

感谢楼主分享

无悔呀

感谢楼主分享教程