由高斯投影坐标反算成经纬度

最新推荐文章于 2026-04-04 09:32:37 发布

转载最新推荐文章于 2026-04-04 09:32:37 发布 · 3.8k 阅读

11 ·

本内容遵循CC 4.0 BY-SA版权协议

原文链接：https://javajgs.com/archives/17789

标签

#线性代数 #几何学 #矩阵

//  由高斯投影坐标反算成经纬度
public static double[] GaussToBL(double X, double Y)
    {
        double[] output = new double[2];
        double longitude1,latitude1, longitude0, X0,Y0, xval,yval;
        //NN曲率半径，测量学里面用N表示
        //M为子午线弧长，测量学里用大X表示 
        //fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示
        //R为底点所对的曲率半径，测量学里用Nf表示
        double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
        iPI = 0.0174532925199433; 3.1415926535898/180.0;
        a=6378137.0; f=1/298.257222101; //CGCS2000坐标系参数
        //a=6378137.0; f=1/298.2572236; //wgs84坐标系参数
        longitude0 = 114.0;//中央子午线
        longitude0 = longitude0 * iPI ; //中央经线

        X0 = 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+10*C-4*C*C-9*ee)*D*D*D*D/24
                +(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);
        //转换为度 DD
        output[0] = longitude1 / iPI;
        output[1] = latitude1 / iPI;
        return output;
    }