高斯坐标正反算示例代码

高斯坐标正反算是地理信息系统（GIS）和大地测量中常用的坐标转换方法。高斯投影是一种横轴墨卡托投影，常用于大比例尺地图的绘制。高斯坐标正算是将地理坐标（经纬度）转换为平面直角坐标（高斯坐标），而高斯坐标反算则是将平面直角坐标转换回地理坐标。

下面是一个简单的Python示例代码，展示了高斯坐标的正反算过程：

import math

# 常量定义
a = 6378137.0  # 地球长半轴
b = 6356752.314245  # 地球短半轴
e_sq = 1 - (b ** 2) / (a ** 2)  # 第一偏心率平方
k0 = 0.9996  # 比例因子
longitude_origin = 0  # 中央子午线

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

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

def geodetic_to_geocentric(latitude, longitude):
    # 将地理坐标转换为地心坐标
    latitude_rad = rad(latitude)
    longitude_rad = rad(longitude)
    
    N = a / math.sqrt(1 - e_sq * math.sin(latitude_rad) ** 2)
    
    X = (N + 0) * math.cos(latitude_rad) * math.cos(longitude_rad)
    Y = (N + 0) * math.cos(latitude_rad) * math.sin(longitude_rad)
    Z = (N * (1 - e_sq)) * math.sin(latitude_rad)
    
    return X, Y, Z

def geocentric_to_geodetic(X, Y, Z):
    # 将地心坐标转换为地理坐标
    p = math.sqrt(X ** 2 + Y ** 2)
    latitude_rad = math.atan2(Z, p * (1 - e_sq))
    
    for _ in range(5):
        N = a / math.sqrt(1 - e_sq * math.sin(latitude_rad) ** 2)
        height = p / math.cos(latitude_rad) - N
        latitude_rad = math.atan2(Z, p * (1 - e_sq * (N / (N + height))))
    
    latitude = deg(latitude_rad)
    longitude = deg(math.atan2(Y, X))
    
    return latitude, longitude

def geodetic_to_gauss(latitude, longitude):
    # 将地理坐标转换为高斯坐标
    latitude_rad = rad(latitude)
    longitude_rad = rad(longitude)
    
    e_prime_sq = e_sq / (1 - e_sq)
    
    N = a / math.sqrt(1 - e_sq * math.sin(latitude_rad) ** 2)
    T = math.tan(latitude_rad) ** 2
    C = e_prime_sq * math.cos(latitude_rad) ** 2
    A = (longitude_rad - rad(longitude_origin)) * math.cos(latitude_rad)
    
    M = a * ((1 - e_sq / 4 - 3 * e_sq ** 2 / 64 - 5 * e_sq ** 3 / 256) * latitude_rad - (3 * e_sq / 8 + 3 * e_sq ** 2 / 32 + 45 * e_sq ** 3 / 1024) * math.sin(2 * latitude_rad) + (15 * e_sq ** 2 / 256 + 45 * e_sq ** 3 / 1024) * math.sin(4 * latitude_rad) - (35 * e_sq ** 3 / 3072) * math.sin(6 * latitude_rad))
    
    easting = k0 * N * (A + (1 - T + C) * A ** 3 / 6 + (5 - 18 * T + T ** 2 + 72 * C - 58 * e_prime_sq) * A ** 5 / 120)
    northing = k0 * (M + N * math.tan(latitude_rad) * (A ** 2 / 2 + (5 - T + 9 * C + 4 * C ** 2) * A ** 4 / 24 + (61 - 58 * T + T ** 2 + 600 * C - 330 * e_prime_sq) * A ** 6 / 720))
    
    return easting, northing

def gauss_to_geodetic(easting, northing):
    # 将高斯坐标转换为地理坐标
    e1 = (1 - math.sqrt(1 - e_sq)) / (1 + math.sqrt(1 - e_sq))
    
    M = northing / k0
    mu = M / (a * (1 - e_sq / 4 - 3 * e_sq ** 2 / 64 - 5 * e_sq ** 3 / 256))
    
    e1 = (1 - math.sqrt(1 - e_sq)) / (1 + math.sqrt(1 - e_sq))
    J1 = (3 * e1 / 2 - 27 * e1 ** 3 / 32)
    J2 = (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32)
    J3 = (151 * e1 ** 3 / 96)
    J4 = (1097 * e1 ** 4 / 512)
    
    latitude_rad = mu + J1 * math.sin(2 * mu) + J2 * math.sin(4 * mu) + J3 * math.sin(6 * mu) + J4 * math.sin(8 * mu)
    
    C1 = e_prime_sq * math.cos(latitude_rad) ** 2
    T1 = math.tan(latitude_rad) ** 2
    N1 = a / math.sqrt(1 - e_sq * math.sin(latitude_rad) ** 2)
    R1 = a * (1 - e_sq) / math.pow(1 - e_sq * math.sin(latitude_rad) ** 2, 1.5)
    D = easting / (N1 * k0)
    
    latitude = latitude_rad - (N1 * math.tan(latitude_rad) / R1) * (D ** 2 / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 ** 2 - 9 * e_prime_sq) * D ** 4 / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 ** 2 - 252 * e_prime_sq - 3 * C1 ** 2) * D ** 6 / 720)
    longitude = rad(longitude_origin) + (D - (1 + 2 * T1 + C1) * D ** 3 / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 ** 2 + 8 * e_prime_sq + 24 * T1 ** 2) * D ** 5 / 120) / math.cos(latitude_rad)
    
    latitude = deg(latitude)
    longitude = deg(longitude)
    
    return latitude, longitude

# 示例使用
latitude = 39.909736  # 纬度
longitude = 116.413772  # 经度

easting, northing = geodetic_to_gauss(latitude, longitude)
print(f"高斯坐标: {easting}, {northing}")

latitude_new, longitude_new = gauss_to_geodetic(easting, northing)
print(f"地理坐标: {latitude_new}, {longitude_new}")