高斯投影反算

爱吃辣条的小王

已于 2024-06-11 17:16:21 修改

阅读量274

点赞数 1

文章标签： python

于 2024-06-11 16:44:05 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_47396624/article/details/139602918

版权

# coding=utf-8
import math

"""2000国家大地坐标系（cgcs2000）椭球常数"""  
a = 6378137.0  # 椭球长半轴
b = 6356752.31414036  # 椭球短半轴
f = 1.0 / 298.257222101  # 椭球扁率
e1 = 0.0066943800229  # 椭球第一偏心率 e^2
ee = 0.00673949677548  # 椭球第二偏心率 e`^2
c = 6399593.62586  # 极点子午圈曲率半径

def xytoBl():
    m0 = a * (1 - e1)
    m2 = 1.5 * e1 * m0
    m4 = 5.0 / 4.0 * e1 * m2
    m6 = 7.0 / 6.0 * e1 * m4
    m8 = 9.0 / 8.0 * e1 * m6

    a0 = m0 + 0.5 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8
    a2 = 0.5 * m2 + 0.5 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8
    a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8
    a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8
    a8 = 1.0 / 128.0 * m8

    X = float(input("请输入x坐标值："))
    Y = float(input("请输入y坐标值：")) - 36000000 - 500000
    l0 = float(input("请输入中央子午线：")) / 180.0 * math.pi

# 迭代求取低点纬度
    B = [0] * 20
    f = [0] * 20
    B[0] = float(X) / a0
    f[0] = -(a2 / 2.0 * math.sin(2 * B[0]) + a4 / 4.0 * math.sin(4.0 * B[0]) - a6 / 6.0 * math.sin(
        6.0 * B[0]) + a8 / 8.0 * math.sin(8.0 * B[0]))
    deltad = B[1] - B[0]
    for i in range(10):
        f[i] = -(a2 / 2.0 * math.sin(2 * B[i]) + a4 / 4.0 * math.sin(4.0 * B[i]) - a6 / 6.0 * math.sin(
            6.0 * B[i]) + a8 / 8.0 * math.sin(8.0 * B[i]))
        B[i + 1] = (X - f[i]) / a0
        deltad = B[i + 1] - B[i]
        if deltad < pow(10, -10):
            i += 1
            break

    """参数"""
    tf = math.tan(B[i])
    mf = c / pow(math.sqrt(1 + ee * pow(math.cos(B[i]), 2)), 3)
    nf = c / math.sqrt(1 + ee * pow(math.cos(B[i]), 2))
    n = math.sqrt(ee) * math.cos(B[i])

    # longitude纬度计算公式
    latitude = B[i] - (tf / (2.0 * mf) * Y * (Y / nf)) * (1.0 - (1.0 / 12.0) * (
            5.0 + 3.0 * pow(tf, 2) + pow(n, 2) - 9.0 * pow(n, 2) * pow(tf, 2)) *
                                                          pow(Y / nf, 2) + 1.0 / 360.0 * (
                                                                  61.0 + 90.0 * pow(tf, 2) + 45.0 * pow(tf, 4)) *
                                                          pow(Y / nf, 4))
    # longitude 纬度计算公式
    longitude = 1.0 / math.cos(B[i]) * (Y / nf) * (
            1.0 - 1.0 / 6.0 * (1.0 + 2.0 * pow(tf, 2) + pow(n, 2)) * pow(Y / nf, 2) + 1.0 / 120.0 * (
            5.0 + 28.0 * pow(tf, 2) + 24.0 * pow(tf, 4) + 6.0 * pow(n, 2) + 8.0 * pow(n, 2) * pow(tf, 2)) * pow(Y / nf,
                                                                                                                4))

    # 将经纬度的弧度转换为度
    lat = latitude * 180 / 3.14159265358979
    long = (longitude + l0) * 180 / 3.14159265358979

    return lat, long

2000坐标系高斯反算；

计算过程中的角度均为弧度；

最后将弧度转换为角度，角度 = 弧度 * 180 / pi