大地坐标与地心坐标相互转换 （WGS84，西安80，北京54， China200）C++

Date 2019/04/09 Add by WJB

我们先了解一下大地坐标系。

地心坐标：地心坐标系（geocentric coordinate system ）以地球质心为原点建立的空间直角坐标系，或以球心与地球质心重合的地球椭球面为基准面所建立的大地坐标系；以地球质心(总椭球的几何中心)为原点的大地坐标系。通常分为地心空间直角坐标系(以x，y，z为其坐标元素)。

大地坐标系：大地坐标（Geodetic coordinate）是大地测量中以参考椭球面为基准面的坐标，地面点P的位置用大地经度L、大地纬度B和大地高H表示。大地坐标多应用于大地测量学，测绘学等。

参考椭球体是一个数学上定义的地球表面，它近似于大地水准面。 由于其相对简单，参考椭球是大地控制网计算和显示点坐标（如纬度，经度和海拔）的首选的地球表面的几何模型。通常所说地球的形状和大小，实际上就是以参考椭球体的长半轴、短半轴和扁率来表示的。

常见椭球体：

北京54：克拉索夫椭球体；西安80：1975国际椭球体，China2000：2000中国大地坐标系。

坐标转换公式：

 

代码：地心转大地

#include<math.h>

double targetH, targetB, targetL;
 void XYZtoBLH(double X, double Y, double Z, double aAxis, double bAxis) {
	double e1 = (pow(aAxis, 2) - pow(bAxis, 2)) / pow(aAxis, 2);
	double e2 = (pow(aAxis, 2) - pow(bAxis, 2)) / pow(bAxis, 2);

	double S = sqrt(pow(X, 2) + pow(Y, 2));
	double cosL = X / S;
	double B = 0;
	double L = 0;

	L = acos(cosL);
	L =fabs(L);

	double tanB = Z / S;
	B = atan(tanB);
	double c = aAxis * aAxis / bAxis;
	double preB0 = 0.0;
	double ll = 0.0;
	double N = 0.0;
	//迭代计算纬度
	do {
		preB0 = B;
		ll =pow(cos(B), 2) * e2;
		N = c / sqrt(1 + ll);

		tanB = (Z + N * e1 *sin(B)) / S;
		B = atan(tanB);
	} while (fabs(preB0 - B) >= 0.0000000001);

	ll = pow(cos(B), 2) * e2;
	N = c /sqrt(1 + ll);

	targetH = s /cos(B) - N * (1 - e1);
	targetB = B * 180 / M_PI;
	targetL = L * 180 / M_PI;
}

代码：地心坐标转大地坐标

 double targetX, targetY, targetZ;
 void BLHtoXYZ(double B, double L, double H, double aAxis, double bAxis) {
	 double dblD2R = M_PI / 180;
	 double e1 = sqrt(pow(aAxis, 2) - pow(bAxis, 2)) / aAxis;

	 double N = aAxis / sqrt(1.0 - pow(e1, 2) * pow(sin(B * dblD2R), 2));
	 targetX = (N + H) * cos(B * dblD2R) * cos(L * dblD2R);
	 targetY = (N + H) * cos(B * dblD2R) * sin(L * dblD2R);
	 targetZ = (N * (1.0 - pow(e1, 2)) + H) * sin(B * dblD2R);
 }