# 高斯坐标正反算算法 # 在正算的时候，y的第一个数字表示是用几度带换算的，第二个数字和第三个数字表示带号 # 该程序是根据武汉大学出版社出版的第三版的《大地测量学基础》编写 # 用于课程实验学习


下面以克氏椭球参数3度带为例进行编程，需要详细的所有的3度带和6度带的各种椭球参数的高斯投影正反算程序代码，可以在本人的上传资料中下载

正算代码：采用的是克氏椭球参数

import math as m

P = 180 / m.pi * 3600  # 1弧度的秒数常数值，1弧度=180/pi=57.29578*3600=206264.80624709636

def GetGS_ZS_KS_3(b, L):
    # 克氏椭球高斯正算
    l = (L - Get3DL0(L)[0]) * 3600 / P
    cos2b = m.pow(m.cos(m.radians(b)), 2)
    N = 6399698.902 - (21562.267 - (108.973 - 0.612 * cos2b) * cos2b) * cos2b  # P.187
    a0 = 32140.404 - (135.3302 - (0.7092 - 0.0040 * cos2b) * cos2b) * cos2b
    a4 = (0.25 + 0.00252 * cos2b) * cos2b - 0.04166
    a6 = (0.166 * cos2b - 0.084) * cos2b
    a3 = (0.3333333 + 0.001123 * cos2b) * cos2b - 0.1666667
    a5 = 0.0083 - (0.1667 - (0.1968 + 0.0040 * cos2b) * cos2b) * cos2b

    x = 6367558.4969 * b * 3600 / P - (a0 - (0.5 + (a4 + a6 * m.pow(l, 2)) * m.pow(l, 2)) * m.pow(l, 2) * N) * m.sin(
        m.radians(b)) * m.cos(m.radians(b))
    y = (1 + (a3 + a5 * m.pow(l, 2)) * m.pow(l, 2)) * l * N * m.cos(m.radians(b)) + Get3DL0(L)[
        1] * 1000000 + 500000 + 300000000
    # 坐标加3度带带号，y坐标西移500000M,加3000000000在开头第一个数字为3，表示为是三度带的转换

    return x, y

反算代码：采用的是克氏椭球参数

def GetGS_FS_KS_3(x, y, L0):
    # 克氏椭球高斯反算
    y = y - L0 / 3 * 1000000 - 500000 - 300000000  # 去掉带号和y坐标恢复东移500000M
    bT = x / 6367558.4969  # bT=x/6367558.4969*P 得到是秒数，不乘P直接得到弧度数值
    beta = x / 6367558.4969 * P  # 得到是秒数
    cos2bT = m.pow(m.cos(bT), 2)
    Bf = (beta + (50221746 + (293622 + (2350 + 22 * cos2bT) * cos2bT) * cos2bT) * m.pow(10, -10) * m.cos(bT) * m.sin(bT)*P)/P
    cos2Bf = m.pow(m.cos(Bf), 2)
    Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * cos2Bf) * cos2Bf) * cos2Bf
    Z = y / (Nf * m.cos(Bf))
    b2 = (0.5 + 0.003369 * cos2Bf) * m.sin(Bf) * m.cos(Bf)
    b3 = 0.333333 - (0.166667 - 0.001123 * cos2Bf) * cos2Bf
    b4 = 0.25 + (0.16161 + 0.00562 * cos2Bf) * cos2Bf
    b5 = 0.2 - (0.1667 - 0.0088 * cos2Bf) * cos2Bf

    b = (Bf * P - (1 - (b4 - 0.12 * m.pow(Z, 2)) * m.pow(Z, 2)) * m.pow(Z, 2) * b2 * P) / 3600.0
    l = ((1 - (b3 - b5 * m.pow(Z, 2)) * m.pow(Z, 2)) * Z * P) / 3600.0
    L = L0 + l

    return b, L

测试数据如下：

正算：

print("程序运行示例如下：")
print("高斯投影正算，克氏椭球参数，用3度带表示，最前面有一个3表示带的种类(小数点前有九位数)：")
print("输入B,L")
print(GetGS_ZS_KS_3(30.5, 114 + 1 / 3))

反算：

print("高斯投影反算，克氏椭球参数，带有三度带带号，最前面还有一个3表示带的种类(小数点前有九位数)：")
print("第一个和第二个参数分别是x,y，第三个输入的参数是中央子午线")
print(GetGS_FS_KS_3(3375648.958456849, 338532000.2685031, 114))

运行结果如下：

正算：

程序运行示例如下：
高斯投影正算，克氏椭球参数，用3度带表示，最前面有一个3表示带的种类(小数点前有九位数)：
输入B,L
(3375648.958456849, 338532000.2685031)

反算：

高斯投影反算，克氏椭球参数，带有三度带带号，最前面还有一个3表示带的种类(小数点前有九位数)：
第一个和第二个参数分别是x,y，第三个输入的参数是中央子午线
(30.500000005448364, 114.33333333335709)