北京54或国家80或CGCS2000转WGS84坐标系的程序实现方法

介绍

在工程测量中，经常使用的坐标系如北京54坐标系、国家80(也叫西安80)坐标系、CGCS2000坐标系。而在终端设备或定位平台中，一般使用WGS84坐标系。这样的情况下，就需要一种北京54、国家80、CGCS2000转WGS84坐标系的通用程序方法。
而我们使用的就是布尔莎七参数法来进行的程序转换。通过使用C++程序，从而得到平台需要的wgs84坐标。

布尔莎七参数

标准的七参数转换方法，使用X、Y、Z平移，X、Y、Z旋转，K尺度，这七个参数，在已知三个以上点的情况下，来计算七参数模型。
具体七参数计算过程，可参考COORD GM软件。

c++、java程序实现坐标转换的过程

c++、java的转换过程，主要依托的参数为七参数，以及原椭球坐标系以及目标椭球坐标系的长半轴、短半轴、扁率、中央子午线。在获取到原始坐标后，经过一系列的坐标转换算法，即可得到经北京54、国家80、CGCS2000坐标系到WGS84坐标系的坐标转换值。

C++代码示例

Point Trans::gauss_geodetic2plane(Point point, double a, double b,double f, int ZoneWide, int Lo)
{
	//计算中央经线和带号
	int ProjNo = 0;
	if (Lo == DEFFAULT_LONTITUDE)//未输入中央经线
	{
		ProjNo = ((int)point.x) / 6 + 1;//带号
		Lo = (ProjNo)* ZoneWide - ZoneWide / 2;//中央经线
	}
	else{
		ProjNo = (Lo + ZoneWide / 2) / ZoneWide;
	}
	double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
	double  e2, ee, NN, T, C, A, M, iPI;
	
	longitude0 = toRadians(Lo);
	longitude1 = toRadians(point.x);
	latitude1 = toRadians(point.y);
	
	e2 = 2 * f - f * f;
	ee = e2 * (1.0 - e2);
	NN = a / sqrt(1.0 - e2 * sin(latitude1) * sin(latitude1));
	T = tan(latitude1) * tan(latitude1);
	C = ee * cos(latitude1) * cos(latitude1);
	A = (longitude1 - longitude0) * 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) * sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * 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 * 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);

	Point pointResult;
	pointResult.x = xval + 1000000 * ProjNo + 500000;
	pointResult.y = yval;
	pointResult.z = point.z;
	return pointResult;
}

Java代码示例

    /**
     * 高斯投影 大地坐标->平面坐标
     */
    Point gauss_geodetic2plane(Point point, double a, double b,double f, int ZoneWide, int Lo)
    {
        //计算中央经线和带号
        int ProjNo = 0;
        if (Lo == DEFFAULT_LONTITUDE)//未输入中央经线
        {
            ProjNo = ((int)point.x) / 6 + 1;//带号
            Lo = (ProjNo)* ZoneWide - ZoneWide / 2;//中央经线
        }
        else{
            ProjNo = (Lo + ZoneWide / 2) / ZoneWide;
        }
        double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
        double  e2, ee, NN, T, C, A, M, iPI;

        longitude0 = toRadians(Lo);
        longitude1 = toRadians(point.x);
        latitude1 = toRadians(point.y);

        e2 = 2 * f - f * f;
        ee = e2 * (1.0 - e2);
        NN = a / sqrt(1.0 - e2 * sin(latitude1) * sin(latitude1));
        T = tan(latitude1) * tan(latitude1);
        C = ee * cos(latitude1) * cos(latitude1);
        A = (longitude1 - longitude0) * 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) * sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * 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 * 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);

        Point pointResult = new Point();
        pointResult.x = xval + 1000000 * ProjNo + 500000;
        pointResult.y = yval;
        pointResult.z = point.z;
        return pointResult;

    }

扩展

实际上，该方法同样支持北京54、国家80、CGCS2000、WGS84坐标系之间的互相转换，该方法是通用的。

如需相关Java及C++源码，请联系我