高斯投影正反算Python源码

# 高斯坐标正反算算法

import math

P=180/math.pi*3600  #1弧度的秒数常数值，1弧度=180/pi=57.29578*3600=206264.80624709636

def GetGSZS54（BL）: # 北京54克拉索夫斯基椭球参数  BL 为十进制度，非度分秒格式，python中math.sin 和.cos函数仅支持弧度，需要用math.radians将度数转换为弧度。
    #北京54椭球高斯正算
    l=（L-Get3DL0（L）[0]）*3600/P
    COS2B=math.pow（math.cos（math.radians（B））2）
    N=6399698.902-（21562.267-（108.973-0.612*COS2B）*COS2B）*COS2B
    A0=32140.404-（135.3302-（0.7092-0.0040*COS2B）*COS2B）*COS2B
    A4=（0.25+0.00252*COS2B）*COS2B-0.04166
    A6=（0.166*COS2B-0.084）*COS2B
    A3=（0.3333333+0.001123*COS2B）*COS2B-0.1666667
    A5=0.0083-（0.1667-（0.1968+0.0040*COS2B）*COS2B）*COS2B

    x=6367558.4969*B*3600/P-（A0-（0.5+（A4+A6*math.pow（l2））*math.pow（l2））*math.pow（l2）*N）*math.sin（math.radians（B））*math.cos（math.radians（B））
    y=（1+（A3+A5*math.pow（l2））*math.pow（l2））*l*N*math.cos（math.radians（B））+Get3DL0（L）[1]*1000000+500000  #坐标加3度带带号，y坐标西移500KM

    return（xy）


def GetGSZS80（BL）:     # 西安80  国际1975国际椭球参数
    #西安80国际1975椭球高斯正算
    l=（L-Get3DL0（L）[0]）*3600/P
    COS2B=math.pow（math.cos（math.radians（B））2）
    N=6399596.652-（21565.045-（108.996-0.603*COS2B）*COS2B）*COS2B
    A0=32144.5189-（135.3646-（0.7034-0.0041*COS2B）*COS2B）*COS2B
    A4=（0.25+0.00253*COS2B）*COS2B-0.04167
    A6=（0.167*COS2B-0.083）*COS2B
    A3=（0.3333333+0.001123*COS2B）*COS2B-0.1666667
    A5=0.00878-（0.1702-0.20382*COS2B）*COS2B

    x=6367452.1328*B*3600/P-（A0-（0.5+（A4+A6*math.pow（l2））*math.pow（l2））*math.pow（l2）*N）*math.cos（math.radians（B））*math.sin（math.radians（B））
    y=（1+（A3+A5*math.pow（l2））*math.pow（l2））*l*N*math.cos（math.radians（B））+Get3DL0（L）[1]*1000000+500000  #坐标加3度带带号，y坐标西移500KM


    return（xy）

def GetGSFS54（xyL0）:
     #北京54椭球高斯反算
    y=y-L0/3*1000000-500000  #去掉带号和y坐标恢复东移500KM
    BT=x/6367558.4969  #  BT=x/6367558.4969*P 得到是秒数，不乘P直接得到弧度数值
    bt=x/6367558.496