Python实现高斯投影正反算

利用Python实现不同椭球坐标系高斯坐标与地理坐标相互转换。

#由经纬度反算成高斯投影坐标（北京54坐标系）
#作者：小熊（初学Python的菜鸟\（￣▽￣）／）

import math

def GTXY1（）:

    def GaussToXY（longitudelatitude）:

        iPI = 3.1415926535898/180.0
        ZoneWide = 6; #6度带宽

        #54年北京坐标系参数

        a=6378245.0

        f=1.0/298.3

        #80年西安坐标系参数
        #a=6378140.0;
        #f=1/298.257;


        ProjNo = （int）（longitude / ZoneWide）

        longitude0 = ProjNo * ZoneWide + ZoneWide / 2 #中央子午线计算
        longitude0 = longitude0 * iPI ;
        latitude0=0;

        longitude1 = longitude * iPI #经度转换为弧度

        latitude1 = latitude * iPI #纬度转换为弧度

        e2=2*f-f*f #第一偏心率

        ee=e2*（1.0-e2）

        NN=a/math.sqrt（1.0-e2*math.sin（latitude1）*math.sin（latitude1））

        T=math.tan（latitude1）*math.tan（latitude1）

        C=ee*math.cos（latitude1）*math.cos（latitude1）

        A=（longitude1-longitude0）*math.cos（latitude1）

        M=a*（（1-e2/4-3*e2*e2/64-5*e2*e2*e2/256）*latitude1-（3*e2/8+3*e2*e2/32+45*e2*e2*e2/1024）*math.sin（2*latitude1）+（15*e2*e2/256+45*e2*e2*e2/1024）*math.sin（4*latitude1）-（35*e2*e2*e2/3072）*math.sin（6*latitude1））

        xval = NN*（A+（1-T+C）*A*A*A/6+（5-18*T+T*T+14*C-58*ee）*A*A*A*A*A/120）

        yval = M+NN*math.tan（latitude1）*（A*A/2+（5-T+9*C+4*C*C）*A*A*A*A/24+（61-58*T+T*T+270*C-330*ee）*A*A*A*A*A*A/720）

        X0 = 1000000L*（ProjNo+1）+500000L 

        Y0 = 0

        xval = xval+X0

        yval = yval+Y0

        return xval yval

    longitude = input（‘输入经度:‘）

    print ‘‘

    latitude = input（‘输入纬度:‘）

    print ‘‘

    print GaussToXY（longitudelatitude）

def GTXY2（）:

    def XYToGauss（XY）:

        iPI = 3.1415926535898/180.0

        #54年北京坐标系参数
        a = 6378245.0
        f = 1.0/298.3

        #80年西安坐标系参数 
        #a=6378140.0;
        #f=1/298.257;

        ZoneWide = 6

        ProjNo = （int）（X/1000000L） #查找带号

        longitude0 = （ProjNo-1） * ZoneWide + ZoneWide / 2;

        longitude0 = longitude0 * iPI #中央经线

        X0 = ProjNo*1000000L+500000L

        Y0 = 0

        xval = X-X0

        yval = Y-Y0

        e2 = 2*f-f*f

        e1 = （1.0-math.sqrt（1-e2））/（1.0+math.sqrt（1-e2））

        ee = e2/（1-e2）

        M = yval

        u = M/（a*（1-e2/4-3*e2*e2/64-5*e2*e2*e2/256））

        fai = u+（3*e1/2-27*e1*e1*e1/32）*math.sin（2*u）+（21*e1*e1/16-55*e1*e1*e1*e1/32）*math.sin（4*u）+（151*e1*e1*e1/96）*math.sin（6*u）+（1097*e1*e1*e1*e1/512）*math.sin（8*u）

        C = ee*math.cos（fai）*math.cos（fai）

        T = math.tan（fai）*math.tan（fai）

        NN = a/math.sqrt（1.0-e2*math.sin（fai）*math.sin（fai））

        R = a*（1-e2）/math.sqrt（（1-e2*math.sin（fai）*math.sin（fai））*（1-e2*math.sin（fai）*math.sin（fai））*（1-e2*math.sin（fai）*math.sin（fai）））

        D = xval/NN

        #计算经度（Longitude） 纬度（Latitude）

        longitude1 = longitude0+（D-（1+2*T+C）*D*D*D/6+（5-2*C+28*T-3*C*C+8*ee+24*T*T）*D*D*D*D*D/120）/math.cos（fai）;

        latitude1 = fai -（NN*math.tan（fai）/R）*（D*D/2-（5+3*T+

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3614  2017-12-06 16:35  高斯正反算.py