北京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坐标系之间的互相转换,该方法是通用的。