高斯投影坐标与经纬度转换的java源码

最新推荐文章于 2024-08-19 11:41:14 发布

Gis051

最新推荐文章于 2024-08-19 11:41:14 发布

阅读量5k

点赞数 1

分类专栏： gis

gis 专栏收录该内容

4 篇文章 0 订阅

订阅专栏

先收藏再说，转 http://blog.sina.com.cn/s/blog_5656b0020100chbv.html

坐标系术语

大地坐标系：地心坐标系(WGS84)、参心坐标系(北京54，西安80)；

平面坐标系：在大地坐标系的基础上，根据不同投影方式生成的平面坐标系，每种大地坐标系都有对应的平面坐标系；

参考椭球：有WGS84 北京54 西安80大地坐标系，是用经纬度表示的;也有WGS84 北京54 西安80平面坐标系,使用xy表示的；

坐标转换：相同参考椭球大地坐标与平面坐标转换、不同椭球体间的坐标转换(7参数)

球面坐标与平面坐标互转(相同椭球体)

1、两点之间球面距离：勾股定理+余弦定理+弧长的公式导出 (涉及角度使用经纬度转弧度计算)

计算球面上任意两点之间的距离的公式通常有两种：Great-circle distance和Haversine formula，而目前大多数公司都是用的是Haversine公式；原因：Great-circle distance公式用到了大量余弦函数，而两点间距离很短时（比如地球表面上相距几百米的两点），余弦函数会得出0.999…的结果，会导致较大的舍入误差。而Haversine公式采用了正弦函数，即使距离很小，也能保持足够的有效数字。

Haversine公式：半角公式；

Great-circle distance：大圆距离公式

2、两点之间平面距离：球面坐标转为高斯平面坐标，AB²=(x1-x2)²+(y1-y2)²；

//球面(经纬度、大地)坐标 与 平面坐标  互转
//1:100万图幅

public class GaussXYDeal {

    public double[] GaussToBL(double X, double Y)   //, double *longitude, double *latitude)
    {
        int ProjNo; int ZoneWide;           带宽
        double[] output = new double[2];
        double longitude1,latitude1, longitude0, X0,Y0, xval,yval;//latitude0,
        double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
        iPI = 0.0174532925199433;           3.1415926535898/180.0;
        //a = 6378245.0; f = 1.0/298.3;     //54年北京坐标系参数
        a=6378140.0; f=1/298.257;           //80年西安坐标系参数
        ZoneWide = 6;                       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+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;
    }

      由经纬度反算成高斯投影坐标
    public double[] BLToGauss(double longitude, double latitude)
    {
        int ProjNo=0; int ZoneWide;             带宽
        double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
        double a,f, e2,ee, NN, T,C,A, M, iPI;
        iPI = 0.0174532925199433;               3.1415926535898/180.0;
        ZoneWide = 6;                           6度带宽
        a=6378245.0; f=1.0/298.3;               //54年北京坐标系参数
        a=6378140.0; f=1/298.257; //80年西安坐标系参数
        ProjNo = (int)(longitude / ZoneWide) ;
        longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
        longitude0 = longitude0 * iPI ;
        latitude0 = 0;
        System.out.println(latitude0);
        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+72*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+600*C-330*ee)*A*A*A*A*A*A/720);
        X0 = 1000000L*(ProjNo+1)+500000L;
        Y0 = 0;

        double[] lonlat = new double[2];
        lonlat[0] = xval+X0;
        lonlat[1] = yval+Y0;
        return lonlat;
    }

     变更后的由经纬度反算成高斯投影坐标, meridian 为mesh中央经度
    public double[] GaussToBLToGauss(double longitude, double latitude, double meridian)
    {
        int ProjNo=0;
        double longitude1,latitude1, X0,Y0, xval,yval;
        double a,f, e2,ee, NN, T,C,A, M, iPI;
        iPI = 0.0174532925199433; 3.1415926535898/180.0;
        a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
        a=6378140.0; f=1/298.257; //80年西安坐标系参数

        meridian = meridian * iPI ;
        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-meridian)*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+72*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+600*C-330*ee)*A*A*A*A*A*A/720);
        X0 = 1000000L*(ProjNo+1)+500000L;
        Y0 = 0;
        double[] lonlat = new double[2];
        lonlat[0] = xval+X0;
        lonlat[1] = yval+Y0;
        return lonlat;
    }

    public static void main(String[] args) throws Exception {
        GaussXYDeal g = new GaussXYDeal();
        double[] xy = g.BLToGauss(113.26846,23.55979);
        System.out.println();

        double[] lonlat = g.GaussToBL(xy[0], xy[1]);
        System.out.println("lon: "+lonlat[0]);
        System.out.println("lat: "+lonlat[1]);

        System.out.println();
        double[] a = g.GaussToBLToGauss(113.26846,23.55979, 113.3125);
        System.out.println("x: "+a[0]);
        System.out.println("y: "+a[1]);
    }

}