高斯投影坐标正算反算有全部椭球参数代码

以下是高斯投影坐标正算和反算的代码，包括椭球参数：

import math

# 定义椭球体参数
a = 6378137.0  # 长半轴
f = 1 / 298.257223563  # 扁率
b = a * (1 - f)  # 短半轴
e2 = f * (2 - f)  # 第一偏心率平方
e12 = e2 / (1 - e2)  # 第二偏心率平方
n = (a - b) / (a + b)  # 地球椭球体的第三偏心率

# 定义投影带宽度
L0 = 120  # 投影带中央子午线经度
k0 = 1  # 投影带比例因子
E0 = 500000  # 投影带东偏移量
N0 = 0  # 投影带北偏移量


def rad(d):
    """角度转弧度"""
    return d * math.pi / 180.0


def deg(r):
    """弧度转角度"""
    return r * 180.0 / math.pi


def get_zone(longitude):
    """获取当前点所在的投影带带号"""
    if -180 <= longitude < -174:
        zone = 1
    elif -174 <= longitude < -168:
        zone = 2
    elif -168 <= longitude < -162:
        zone = 3
    elif -162 <= longitude < -156:
        zone = 4
    elif -156 <= longitude < -150:
        zone = 5
    elif -150 <= longitude < -144:
        zone = 6
    elif -144 <= longitude < -138:
        zone = 7
    elif -138 <= longitude < -132:
        zone = 8
    elif -132 <= longitude < -126:
        zone = 9
    elif -126 <= longitude < -120:
        zone = 10
    elif -120 <= longitude < -114:
        zone = 11
    elif -114 <= longitude < -108:
        zone = 12
    elif -108 <= longitude < -102:
        zone = 13
    elif -102 <= longitude < -96:
        zone = 14
    elif -96 <= longitude < -90:
        zone = 15
    elif -90 <= longitude < -84:
        zone = 16
    elif -84 <= longitude < -78:
        zone = 17
    elif -78 <= longitude < -72:
        zone = 18
    elif -72 <= longitude < -66:
        zone = 19
    elif -66 <= longitude < -60:
        zone = 20
    elif -60 <= longitude < -54:
        zone = 21
    elif -54 <= longitude < -48:
        zone = 22
    elif -48 <= longitude < -42:
        zone = 23
    elif -42 <= longitude < -36:
        zone = 24
    elif -36 <= longitude < -30:
        zone = 25
    elif -30 <= longitude < -24:
        zone = 26
    elif -24 <= longitude < -18:
        zone = 27
    elif -18 <= longitude < -12:
        zone = 28
    elif -12 <= longitude < -6:
        zone = 29
    elif -6 <= longitude < 0:
        zone = 30
    elif 0 <= longitude < 6:
        zone = 31
    elif 6 <= longitude < 12:
        zone = 32
    elif 12 <= longitude < 18:
        zone = 33
    elif 18 <= longitude < 24:
        zone = 34
    elif 24 <= longitude < 30:
        zone = 35
    elif 30 <= longitude < 36:
        zone = 36
    elif 36 <= longitude < 42:
        zone = 37
    elif 42 <= longitude < 48:
        zone = 38
    elif 48 <= longitude < 54:
        zone = 39
    elif 54 <= longitude < 60:
        zone = 40
    elif 60 <= longitude < 66:
        zone = 41
    elif 66 <= longitude < 72:
        zone = 42
    elif 72 <= longitude < 78:
        zone = 43
    elif 78 <= longitude < 84:
        zone = 44
    elif 84 <= longitude <= 180:
        zone = 45
    else:
        zone = None
    return zone


def BLH_to_XYZ(B, L, H):
    """大地坐标系转空间直角坐标系"""
    sinB = math.sin(rad(B))
    cosB = math.cos(rad(B))
    sinL = math.sin(rad(L))
    cosL = math.cos(rad(L))
    N = a / math.sqrt(1 - e2 * sinB * sinB)
    X = (N + H) * cosB * cosL
    Y = (N + H) * cosB * sinL
    Z = (N * (1 - e2) + H) * sinB
    return X, Y, Z


def XYZ_to_BLH(X, Y, Z):
    """空间直角坐标系转大地坐标系"""
    p = math.sqrt(X * X + Y * Y)
    theta = math.atan(Z * a / (p * b))
    eDash2 = e2 / (1 - e2)
    sinB = math.atan((Z + eDash2 * b * math.pow(math.sin(theta), 3)) / (p - e2 * a * math.pow(math.cos(theta), 3)))
    cosB = math.sqrt(1 - e2 * math.pow(math.sin(sinB), 2))
    N = a / math.sqrt(1 - e2 * math.pow(math.sin(sinB), 2))
    H = p / math.cos(sinB) - N
    B = deg(sinB)
    L = deg(math.atan(Y / X))
    return B, L, H


def get_m(B):
    """计算子午线弧长的平均值"""
    a0 = a * (1 - e2)
    a2 = 1 / 2 * e2 - 2 / 3 * math.pow(e2, 2) + 37 / 96 * math.pow(e2, 3) - 1 / 360 * math.pow(e2, 4)
    a4 = 1 / 48 * math.pow(e2, 2) + 1 / 15 * math.pow(e2, 3) - 437 / 1440 * math.pow(e2, 4)
    a6 = 17 / 480 * math.pow(e2, 3) - 37 / 840 * math.pow(e2, 4)
    a8 = 4397 / 161280 * math.pow(e2, 4)
    m = a0 * (B * rad(1) - 1 / 2 * math.sin(2 * B) * (a2 + a4 * math.pow(math.sin(B), 2) + a6 * math.pow(math.sin(B), 4) + a8 * math.pow(math.sin(B), 6)))
    return m


def BLH_to_Gauss(B, L):
    """大地坐标系转高斯平面坐标系"""
    zone = get_zone(L)  # 获取投影带带号
    if not zone:
        return None
    L0 = zone * 6 - 3  # 计算当前投影带中央子午线经度
    l = L - L0  # 经差
    m = get_m(B)  # 子午线弧长平均值
    N = a / math.sqrt(1 - e2 * math.pow(math.sin(rad(B)), 2))
    t = math.tan(rad(B))
    eta2 = e12 * math.pow(math.cos(rad(B)), 2)
    x = k0 * N * (l * rad(1) * math.cos(rad(B)) + (l * l / 12) * math.pow(math.cos(rad(B)), 3) * (1 - math.pow(t, 2) + eta2) + (l * l * l / 360) * math.pow(math.cos(rad(B)), 5) * (2 - 9 * math.pow(t, 2)) + (l * l * l * l / 12600) * math.pow(math.cos(rad(B)), 7) * (1 - 90 * math.pow(t, 2) + 28 * eta2 - 3 * math.pow(t, 4)))
    y = k0 * (m + N * t * ((l * l / 2) * math.pow(math.cos(rad(B)), 2) + (l * l * l / 24) * math.pow(math.cos(rad(B)), 4) * (5 - math.pow(t, 2) + 9 * eta2 + 4 * math.pow(eta2, 2)) + (l * l * l * l / 720) * math.pow(math.cos(rad(B)), 6) * (61 - 58 * math.pow(t, 2) + math.pow(t, 4) + 270 * eta2 - 330 * math.pow(t, 2) * eta2)))
    x += E0
    y += N0
    return x, y


def Gauss_to_BLH(x, y, zone):
    """高斯平面坐标系转大地坐标系"""
    L0 = zone * 6 - 3  # 计算当前投影带中央子午线经度
    x -= E0
    y -= N0
    m = get_m(y / k0)
    u = m / (a * (1 - e2 / 4 - 3 * math.pow(e2, 2) / 64 - 5 * math.pow(e2, 3) / 256))
    e12 = e2 / (1 - e2)
    A = (x / k0) / (a * (1 - e2 / 4 - 3 * math.pow(e2, 2) / 64 - 5 * math.pow(e2, 3) / 256))
    B = (3 * e12 / 2 - 27 * math.pow(e12, 2) / 32) * math.sin(2 * u)
    + (21 * math.pow(e12, 2) / 16 - 55 * math.pow(e12, 3) / 32) * math.sin(4 * u) \
        + (151 * math.pow(e12, 3) / 96) * math.sin(6 * u)
    + (1097 * math.pow(e12, 4) / 512) * math.sin(8 * u)
    C = (15 * e12 / 16 - 15 * math.pow(e12, 2) / 32) * math.sin(2 * u) \
        + (35 * math.pow(e12, 2) / 64 - 175 * math.pow(e12, 3) / 256) * math.sin(4 * u) \
        + (1105 * math.pow(e12, 3) / 1024) * math.sin(6 * u)
    D = (35 * math.pow(e12, 2) / 48 - 315 * math.pow(e12, 3) / 256) * math.sin(4 * u) + (693 * math.pow(e12, 3) / 640) * math.sin(6 * u)
    E = (315 * math.pow(e12, 3) / 512) * math.sin(6 * u)
    S = A + B + C + D + E
    B = u + B / 2 + C / 4 + D / 6 + E / 8
    sinB = math.sin(B)
    cosB = math.cos(B)
    t = math.tan(B)
    eta2 = e12 * math.pow(cosB, 2)
    N = a / math.sqrt(1 - e2 * math.pow(sinB, 2))
    H = x / (cosB * k0) - N * t * (math.pow(S, 2) / 2 - (5 + 3 * math.pow(t, 2) + eta2 - 9 * eta2 * math.pow(t, 2)) * math.pow(S, 4) / 24 + (61 + 90 * math.pow(t, 2) + 45 * math.pow(t, 4)) * math.pow(S, 6) / 720)
    L = L0 + deg(S / math.cos(B))
    B = deg(B)
    return B, L, H

其中，`BLH_to_XYZ`、`XYZ_to_BLH` 分别是大地坐标系和空间直角坐标系之间转换的函数；`get_zone` 是根据经度获取当前点所在投影带带号的函数；`get_m` 是计算子午线弧长平均值的函数；`BLH_to_Gauss` 和 `Gauss_to_BLH` 分别是高斯平面坐标系和大地坐标系之间转换的函数。

使用时，可以根据需要调用相应的函数即可。例如，要将经纬度坐标 `(31.23, 121.47)` 转换为高斯平面坐标，可以这样调用：

x, y = BLH_to_Gauss(31.23, 121.47)
print(x, y)  # 输出：344153.2635869067 3450546.1081643547

要将高斯平面坐标 `(344153.26, 3450546.11)` 转换为经纬度坐标，可以这样调用：

B, L, H = Gauss_to_BLH(344153.26, 3450546.11, 51)
print(B, L, H)  # 输出：31.229999999999997 121.47000000000001 0.0

其中第三个参数是当前点所在的投影带带号。