2000国家大地坐标系高斯正算

# coding=utf-8
import math

"""2000国家大地坐标系（cgcs2000）椭球常数"""  # 度=弧度×180°÷π
a = 6378137.0  # 椭球长半轴
b = 6356752.31414036  # 椭球短半轴
f = 1.0 / 298.257222101  # 椭球扁率
e1 = 0.0066943800229  # 椭球第一偏心率 e^2
ee = 0.00673949677548  # 椭球第二偏心率 e`^2
c = 6399593.62586  # 极点子午圈曲率半径

"""相关常数"""
K0 = 1.57048761144159 * math.pow(10.0, -7.0)
K1 = 5.05250178820567 * math.pow(10.0, -3.0)
K2 = 2.98472900956587 * math.pow(10.0, -5.0)
K3 = 2.41626669230084 * math.pow(10.0, -7.0)
K4 = 2.22241238938534 * math.pow(10.0, -9.0)


def BLtoXY():
    """高斯正算 有带号"""
    """------------------------------------------------"""
    B = input("请输入纬度：")  # 输入的必须是弧度
    L = input("请输入经度：")  # 输入的必须是弧度
    L0 = input("请输入中央子午线：")  # 输入的必须是弧度
    """-------------------------------------------------"""
    m0 = a * (1 - e1)
    m2 = 1.5 * e1 * m0
    m4 = 5.0 / 4.0 * e1 * m2
    m6 = 7.0 / 6.0 * e1 * m4
    m8 = 9.0 / 8.0 * e1 * m6

    a0 = m0 + 0.5 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8
    a2 = 0.5 * m2 + 0.5 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8
    a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8
    a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8
    a8 = 1.0 / 128.0 * m8

    N = a / math.sqrt(1 - e1 * math.sin(B) * math.sin(B))

    X = (a0 * B) - a2 / 2.0 * math.sin(2 * B) + a4 / 4.0 * math.sin(4.0 * B) - a6 / 6.0 * math.sin(
        6 * B) + a8 / 8.0 * math.sin(8.0 * B)

    l = L - L0      # 中央子午线
    t = math.tan(B)
    n = math.sqrt(ee) * math.cos(B)

    # x值计算公式
    x = X + N * t * (0.5 * math.cos(B) * math.cos(B) * l * l +
                     1.0 / 24.0 * pow(math.cos(B), 4) * pow(l, 4) * (5.0 - t * t + 9.0 * n * n + 4.0 * pow(n, 4)) +
                     1.0 / 720.0 * pow(math.cos(B),
                                       6) * pow(l,
                                                6) * (
                             61.0 - 58.0 * t * t + pow(t,
                                                       4)))
    # y值计算公式
    y = (N * (math.cos(B) * l + 1.0 / 6.0 * pow(math.cos(B), 3) * pow(l, 3) * (1.0 - t + n * n) +
              1.0 / 120.0 * pow(math.cos(B), 5) * pow(l, 5) * (
                      5.0 - 18.0 * t * t + pow(t, 4) + 14.0 * n * n - 58.0 * n * n * t * t))) + 500000 + 37000000

    return y, x


if __name__ == '__main__':
    try:
        coord = BLtoXY()
        print (coord)

    except Exception as e:
        print (str(e))

输入的经度、纬度和中央子午线必须是弧度。

度转弧度公式如下：

弧度 = 度 / 180 * π