高斯投影坐标正反算python

高斯投影是一种地图投影方式，涉及到高斯投影坐标的正算和反算。下面是Python实现的高斯投影坐标正反算代码示例：

高斯投影坐标正算：

import math

def gauss_forward(lng,lat):
    a = 6378245.0 # 长半轴
    ee = 0.00669342162296594323 # 椭球偏心率平方
    lng0 = 120.0 # 中央经线
    k0 = 1.0 # 比例因子
    dx = 500000.0 # 偏移量

    # 经度转换为弧度
    rad_lng = math.radians(lng)
    # 纬度转换为弧度
    rad_lat = math.radians(lat)

    # 计算子午线弧长
    a0 = a * (1 - ee)
    a2 = a / 2 + (1 - ee) * a / 4 + (5 - 3 * ee) * a / 96
    a4 = (7 - 11 * ee) * a / 48 + (29 - 30 * ee) * a / 192
    a6 = 811 - 13 * ee
    m = a0 * (rad_lat - a2 * math.sin(2 * rad_lat) + a4 * math.sin(4 * rad_lat) - a6 * math.sin(6 * rad_lat))

    # 计算高斯投影坐标
    n = (a - a0) / (a + a0)
    cos_lat = math.cos(rad_lat)
    cos2_lat = cos_lat * cos_lat
    cos4_lat = cos2_lat * cos2_lat
    cos6_lat = cos4_lat * cos2_lat
    cos8_lat = cos4_lat * cos4_lat
    A = (rad_lng - math.radians(lng0)) * cos_lat
    A2 = A * A
    A4 = A2 * A2
    A6 = A4 * A2
    A8 = A6 * A2
    a0 = a * (1 - n + (5 / 4) * (n ** 2 - n ** 3) + (81 / 64) * (n ** 4 - n ** 5))
    a2 = (3 / 2) * (n - n ** 2 + (7 / 8) * (n ** 3 - n ** 4) + (55 / 64) * n ** 5)
    a4 = (15 / 16) * (n ** 2 - n ** 3 + (3 / 4) * (n ** 4 - n ** 5))
    a6 = (35 / 48) * (n ** 3 - n ** 4 + (11 / 16) * n ** 5)
    a8 = (315 / 512) * (n ** 4 - n ** 5)
    x = dx + k0 * a0 * (A + (1 / 2) * a2 * A2 + (1 / 4) * a4 * A4 + (1 / 6) * a6 * A6 + (1 / 8) * a8 * A8)
    y = k0 * (m + a0 * math.tan(rad_lat) * (A2 / 2 + (1 / 3) * a4 * A4 + (1 / 5) * a6 * A6 + (1 / 7) * a8 * A8))

    return x, y

高斯投影坐标反算：

import math

def gauss_inverse(x,y):
    a = 6378245.0 # 长半轴
    ee = 0.00669342162296594323 # 椭球偏心率平方
    lng0 = 120.0 # 中央经线
    k0 = 1.0 # 比例因子
    dx = 500000.0 # 偏移量

    # 计算反算需要的参数
    e1 = (1 - math.sqrt(1 - ee)) / (1 + math.sqrt(1 - ee))
    e12 = e1 ** 2 / (1 - e1 ** 2)
    m = y / k0
    mu = m / (a * (1 - ee / 4 - 3 * ee ** 2 / 64 - 5 * ee ** 3 / 256))
    e2 = e1 ** 2 / (1 - e1 ** 2)
    c1 = e12 * math.cos(mu) ** 2
    t1 = math.tan(mu) ** 2
    n = a / math.sqrt(1 - ee * math.sin(mu) ** 2)
    r = a * (1 - ee) / math.sqrt((1 - ee * math.sin(mu) ** 2) ** 3)

    # 计算经度和纬度
    lng = lng0 + x / (n * k0)
    lat = mu - t1 * (1 + c1) * (x / (n * k0) ** 2) / 2 + t1 * (5 + 3 * t1 + 10 * c1 - 4 * c1 ** 2 - 9 * e2) * (x / (n * k0) ** 4) / 24 - t1 * (61 + 90 * t1 + 298 * c1 + 45 * t1 ** 2 - 252 * e2 - 3 * c1 ** 2) * (x / (n * k0) ** 6) / 720
    lat = math.degrees(lat)
    lng = math.degrees(lng)

    return lng, lat

其中，`lng`和`lat`分别表示经度和纬度，`x`和`y`分别表示高斯投影坐标系下的横坐标和纵坐标。