武汉大学测绘学院——大地坐标与空间直角坐标的相互转换

采用C++语言，实现读取逗号分隔的文本文件，大地坐标和空间直角坐标的相互转换。

1.空间直角坐标与结构体定义

struct XYZ
	{
		double x;
		double y;
		double z;
	};//笛卡尔坐标
	vector<XYZ> P;

2.大地坐标与结构体定义

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

struct BLH
	{
		double B;
		double L;
		double H;
	};//大地坐标
	vector<BLH> B;

3.数据文件读取

数据文件格式如下：

从左到右分别为编号，x，y，z坐标

逐行读取数据的函数

void readfile(string filename)
	{
		ifstream infile(filename);
		if (!infile.is_open())
		{
			cerr << "open error!" << endl;
		}
		string line;
		XYZ P1;
		while (getline(infile, line))
		{
			stringstream ss(line);
			string token;
			if (getline(ss, token, ','))
			{
				int number=stoi(token);
			}
			if (getline(ss, token, ','))
			{
				P1.x = stod(token);
			}
			if (getline(ss, token, ','))
			{
				P1.y = stod(token);
			}
			if (getline(ss, token, ','))
			{
				P1.z = stod(token);
			}
			P.push_back(P1);//读取逗号分隔的文本文件，并添加到容器内
		}
	}

stoi为字符串转换成整数的函数，stod为字符串转换成双精度浮点数的函数

使用动态数组存储数据，方便添加数据

4.空间直角坐标（笛卡尔坐标）到大地坐标的转换算法

当然，在实际计算中还有很多其它的算法，这里选择以上算法。

这一步的难点在于迭代计算，取初值，计算得到B的初值，再用B的值更新ΔZ的值，不断迭代，直到两次的B的差值小于1e-12（实际1e-10即可满足要求），实现代码如下

首先需要有一个常量定义，包括地球偏心率，长半轴等，然后再计算

double a = 6378137;//长半轴
double f = 1 / 298.257223563;//焦距
double b = a * (1 - f);//短半轴
double c = sqrt(a*a - b * b);
double e2 = (a*a - b * b) / (a*a);//第一偏心率的平方
double ep2 = (a*a - b * b) / (b*b);//第二偏心率的平方
void XYZ2BLH()
	{
		double x, y, z;
		BLH B1;
		for (int i = 0; i <= size(P) - 1; i++)
		{
			x = P[i].x;
			y = P[i].y;
			z = P[i].z;
			B1.L = atan2(y, x);//弧度单位的经度
			double N = 0.0;//卯酉圈半径
			double delta_z = e2 * z;//迭代计算的初始值
			double B_current = atan2((z + delta_z), sqrt(x*x + y * y));
			double B_previous = 0.0;
			while (abs(B_current - B_previous) > 1e-12)
			{
				B_previous = B_current;
				N = a / sqrt(1 - e2 * pow(sin(B_previous), 2));
				delta_z = N * e2*sin(B_previous);
				B_current = atan2((z + delta_z), sqrt(x*x + y * y));
			}
			B1.B = B_current;//取最后一次迭代的值
			B1.H = sqrt(pow(x, 2) + pow(y, 2) + pow((z + delta_z), 2)) - N;//计算高程
			B.push_back(B1);
		}
	}

5.大地坐标到空间直角坐标的转换算法

这一步的转换算法相对简单，直接计算即可

vector<XYZ> Pk;
	void BLH2XYZ()
	{
		for (int i = 0; i <= size(B) - 1; i++)
		{
			BLH B1 = B[i];
			double N = a / sqrt(1 - e2 * pow(sin(B1.B), 2));
			double x = (N + B1.H)*cos(B1.B)*cos(B1.L);
			double y = (N + B1.H)*cos(B1.B)*sin(B1.L);
			double z = (N*(1 - e2) + B1.H)*sin(B1.B);
			XYZ P1;
			P1.x = x; P1.y = y; P1.z = z;
			Pk.push_back(P1);
		}
	}

6.结果检验

由于我们写了两种坐标相互转换的函数，所以我们可以将空间直角坐标转换成大地坐标之后的数据再转换回去，如果得到的结果和原来的相同，说明结果正确。

写一个输出函数来做检验

void showdata()
	{
		cout << fixed << setprecision(4);
		for (int i = 0; i <= size(Pk) - 1; i++)
		{
			cout << i + 1 << "," << Pk[i].x << "," << Pk[i].y << "," << Pk[i].z << endl;
		}
	}

输出结果如下：

和原数据一毛一样，说明计算结果正确。

这次的内容就分享到这里了，有问题欢迎大家一起讨论！

完整源代码及测试用例见：

通用时、儒略日、GPS时转换算法和大地坐标与空间直角坐标转换算法资源-CSDN文库

仅需1积分即可获取