高斯投影坐标正算python

以下是一个使用python进行高斯投影坐标正算的例子：

import math

# 定义一些常量
a = 6378137.0  # 长半轴
f = 1 / 298.257223563  # 扁率
L0 = 105  # 中央经线

# 定义一些函数
def rad(d):
    return d * math.pi / 180.0

def get_e2():
    return f * (2 - f)

def get_ep2():
    return get_e2() / (1 - get_e2())

def get_m(B):
    e2 = get_e2()
    a1 = a * (1 - e2)
    a2 = a / (1 - e2 * math.sin(rad(B)) ** 2) ** 0.5
    return a1 / a2

def get_x(B, L):
    l = rad(L - L0)
    m = get_m(B)
    e2 = get_e2()
    ep2 = get_ep2()
    N = a / (1 - e2 * math.sin(rad(B)) ** 2) ** 0.5
    t = math.tan(rad(B)) ** 2
    cosB = math.cos(rad(B))
    x = N * l * cosB + N * (l ** 3) * (1 - t + 5 * (t ** 2 - t + ep2) * (l ** 2) / 24) * cosB ** 3 / 6 + N * (l ** 5) * (5 - 18 * t + (t ** 2) + 72 * ep2 - 58 * ep2 * t) * cosB ** 5 / 120
    return x + 500000

def get_y(B, L):
    l = rad(L - L0)
    m = get_m(B)
    e2 = get_e2()
    ep2 = get_ep2()
    N = a / (1 - e2 * math.sin(rad(B)) ** 2) ** 0.5
    t = math.tan(rad(B)) ** 2
    cosB = math.cos(rad(B))
    y = N * (l ** 2) * math.tan(rad(B)) / 2 + N * (l ** 4) * (5 - t + 9 * ep2 + 4 * (ep2 ** 2) * (l ** 2)) * math.tan(rad(B)) ** 3 / 24 + N * (l ** 6) * (61 - 58 * t + (t ** 2) + 600 * ep2 - 330 * ep2 * t) * math.tan(rad(B)) ** 5 / 720
    return y + 0

# 输入大地坐标
B = 25.0
L = 102.0

# 计算高斯平面坐标
x = get_x(B, L)
y = get_y(B, L)

# 输出结果
print("高斯平面坐标为：({:.2f}, {:.2f})".format(x, y))