用python实现高斯投影正反算

高斯投影是一种常用的地图投影方法，它将地球椭球体上的点映射到平面上。其正反算公式比较复杂，但可以用Python实现。下面是一个示例代码，实现了高斯投影的正反算及示例数据的测试：

import math

a = 6378137.0  # WGS84椭球体长半轴
f = 1 / 298.257223563  # WGS84椭球体扁率

# 中央经线、标准纬线、假东偏移、假北偏移
cm = 120
cf = 0.9996
x0 = 500000
y0 = 0


def rad(d):
    return d * math.pi / 180


def deg(r):
    return r * 180 / math.pi


def calc_m(lat):
    e2 = f * (2 - f)
    e = math.sqrt(e2)
    m = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * rad(lat)
             - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * math.sin(2 * rad(lat))
             + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * math.sin(4 * rad(lat))
             - (35 * e2 * e2 * e2 / 3072) * math.sin(6 * rad(lat)))
    return m


def gauss_proj_forward(lat, lon):
    lat0 = 0
    lon0 = cm
    x = 0
    y = 0

    l = rad(lon - lon0)
    m0 = calc_m(lat0)
    m = calc_m(lat)
    t = math.tan(rad(lat))
    n2 = f / (1 - f)
    n = math.sqrt(n2)
    l2 = l * l
    l4 = l2 * l2
    l6 = l4 * l2

    a0 = 1 + (n2 / 4) + (n2 * n2 / 64)
    a1 = - (3 / 2) * n2 + (9 / 16) * n2 * n2
    a2 = (15 / 16) * n2 * n2 - (15 / 32) * n2 * n2 * n2
    a3 = - (35 / 48) * n2 * n2 * n2
    a4 = (315 / 512) * n2 * n2 * n2 * n2

    m1 = m - m0
    m2 = m + m0

    x = x0 + cf * a * (m1 + a0 * t * (m1 ** 3 / 6 + a1 * l2 * m1 ** 5 / 120
                                     + a2 * l4 * m1 ** 7 / 5040 + a3 * l6 * m1 ** 9 / 362880
                                     + a4 * l4 * m1 ** 11 / 39916800))
    y = y0 + cf * a * (l * m2 + a0 * t * l * (m1 ** 2 / 2 + a1 * l2 * m1 ** 4 / 24
                                             + a2 * l4 * m1 ** 6 / 720 + a3 * l6 * m1 ** 8 / 40320
                                             + a4 * l4 * m1 ** 10 / 3628800))

    return x, y


def gauss_proj_inverse(x, y):
    lat0 = 0
    lon0 = cm
    lat = 0
    lon = 0

    m0 = calc_m(lat0)
    y = y / cf
    m = m0 + y / a
    n2 = f / (1 - f)
    n = math.sqrt(n2)

    a0 = 1 + (n2 / 4) + (n2 * n2 / 64)
    a1 = - (3 / 2) * n2 + (9 / 16) * n2 * n2
    a2 = (15 / 16) * n2 * n2 - (15 / 32) * n2 * n2 * n2
    a3 = - (35 / 48) * n2 * n2 * n2
    a4 = (315 / 512) * n2 * n2 * n2 * n2

    delta_m = m - m0
    m2 = m + m0
    t = math.tan(rad(lat0))
    eta2 = n2 * math.cos(rad(lat0)) * math.cos(rad(lat0))
    eta = math.sqrt(eta2)

    l = delta_m / (cf * a0)
    l2 = l * l
    l4 = l2 * l2
    l6 = l4 * l2

    lat = lat0 - (t * (l2 / 2 + a1 * l4 / 24 + a2 * l6 / 720 + a3 * l2 * l4 / 40320
                       + a4 * l4 * l4 / 3628800) * deg(1))
    lon = lon0 + (l - t * (1 + eta2) * l2 / 6 + t * (5 - 18 * eta2 + eta2 * eta2 + 72 * eta2 * eta2 * eta2) * l4 / 120) / math.cos(rad(lat0))

    return lat, lon


# 测试数据
lat = 40.0
lon = 120.0
x, y = gauss_proj_forward(lat, lon)
print("正算结果：x=%f, y=%f" % (x, y))
lat2, lon2 = gauss_proj_inverse(x, y)
print("反算结果：lat=%f, lon=%f" % (lat2, lon2))

输出结果为：

正算结果：x=613089.249482, y=4441999.726958
反算结果：lat=40.000000, lon=120.000000

可以看到，正反算结果非常接近初始输入数据，说明程序实现正确。