2021-09-09

C++坐标转换（源代码）

//N=a (1-e^2 sin^2 B)^ (-1/2)
//a——椭球长半径；e——椭球第一偏心率；B——大地纬度。
#include<iostream>
#include <cmath>
#define pi 3.1415926535897
using namespace std;
void XYZtoBLH(double X, double Y, double Z) {
	double a, b, f,e2, N[5];//e2是e的平方
	double B[5], L, H[5];
	L = atan(Y / X);
	int chose = 2;
	cout << "chose GRS80 1 or WGS84 2"; cin >> chose;
	//选择长半轴的参数
	if (chose == 1)
	{
			a = 6378137; //% 长半轴
			f = 1 / 298.257222101;// % 扁率
	}
	else
	{
			a = 6378137; //% 长半轴
			f = 1 / 298.257223563;// % 扁率
	}
	b = (1 - f) * a; //% 短半轴
	e2 = 1 - pow(b, 2) / pow(a, 2);
	//迭代结算的初值
	N[0] = a;
	H[0] = pow(X * X + Y * Y + Z * Z, 0.5) - pow(a * b, 0.5);
	B[0] = atan((Z / pow(X * X + Y * Y, 0.5)) * 
			pow((1 - (e2 * N[0]) / (N[0] + H[0])), -1));
	for (int i = 1; i < 5; i++) {
		N[i] = a / pow((1 - e2 * sin(B[i - 1]) * sin(B[i - 1])), 0.5);
		H[i] = pow(X * X + Y * Y, 0.5) / cos(B[i - 1]) - N[i];
		B[i] = atan((Z / pow(X * X + Y * Y, 0.5)) * pow((1 - e2 * N[i] / (N[i] + H[i])), -1));
		if (i ==4) {
			double ph = H[3] - H[2];
			//弧度* (180 / π） = 度数
			double pb = (B[3] - B[2])*(180/pi)*3600;
			if (ph <= 0.001&&pb<0.00001)break;
		}
	}
	cout << 'B' << B[3] * (180 / pi) << endl;
	cout << 'L' << L * (180 / pi) << endl;
	cout << 'H' << H[3] << endl;
}
void BLHtoXYZ(double B, double L, double H) {
	double a, b, f, e2, N;//e2是e的平方
	double X, Y, Z;
	int chose = 2;
	cout << "chose GRS80 1 or WGS84 2"; cin >> chose;
	if (chose == 1)
	{
		a = 6378137; //% 长半轴
		f = 1 / 298.257222101;// % 扁率
	}
	else
	{
		a = 6378137; //% 长半轴
		f = 1 / 298.257223563;// % 扁率
	}
	b = (1 - f) * a; //% 短半轴
	e2 = 1 - pow(b, 2) / pow(a, 2);
	/*cout << "请输入B，" << "L，" << "H:"; cin >> B >> L >> H;*/
	N = a / pow((1 - e2 * (sin(B * pi / 180) * sin(B * pi / 180))), 0.5);
	X = (N + H) * cos(B * pi / 180) * cos(L * pi / 180);
	Y = (N + H) * cos(B * pi / 180) * sin(L * pi / 180);
	Z = (N * (1 - e2) + H) * sin(B * pi / 180);
	cout << 'X' << X << endl;
	cout << 'Y' << Y << endl;
	cout << 'Z' << Z << endl;
}
int main() 
{
	XYZtoBLH(3194419.6, 3199522.3, 9588999.5);
}

直角转换大地

%   intput：
%           (Lon，Lat，H)------------经纬度和大地高，可以为列向量
%
%           ReferenceFrame----------参考框架
%   output：
%           (X,Y,Z)-----------------空间直角坐标
%
function [X,Y,Z]=blh2xyz(Lon,Lat,H,ReferenceFrame)

if nargin<4
    ReferenceFrame='GRS80';
end
switch ReferenceFrame
    %GRS80椭球参数，根据IERS2010得到
    case 'GRS80'
        a=6378137;%长半轴
        f=1/298.257222101;%扁率
    %WGS84椭球参数
    case 'WGS84'
        a=6378137;%长半轴
        f=1/298.257223563;%扁率
end

Lon=Lon*pi/180;
Lat=Lat*pi/180;
b=(1-f)*a;%短半轴
e2=1-b^2/a^2;%偏心率的平方
N=a./sqrt(1-e2*sin(Lat).^2);%卯酉圈曲率半径

%空间直角坐标
X=(N+H).*cos(Lat).*cos(Lon);
Y=(N+H).*cos(Lat).*sin(Lon);
Z=(N*(1-e2)+H).*sin(Lat);

大地转化直角

%   intput：
%           (Lon，Lat，H)------------经纬度和大地高，可以为列向量
%
%           ReferenceFrame----------参考框架
%   output：
%           (X,Y,Z)-----------------空间直角坐标
%
function [X,Y,Z]=blh2xyz(Lon,Lat,H,ReferenceFrame)

if nargin<4
    ReferenceFrame='GRS80';
end
switch ReferenceFrame
    %GRS80椭球参数，根据IERS2010得到
    case 'GRS80'
        a=6378137;%长半轴
        f=1/298.257222101;%扁率
    %WGS84椭球参数
    case 'WGS84'
        a=6378137;%长半轴
        f=1/298.257223563;%扁率
end

Lon=Lon*pi/180;
Lat=Lat*pi/180;
b=(1-f)*a;%短半轴
e2=1-b^2/a^2;%偏心率的平方
N=a./sqrt(1-e2*sin(Lat).^2);%卯酉圈曲率半径

%空间直角坐标
X=(N+H).*cos(Lat).*cos(Lon);
Y=(N+H).*cos(Lat).*sin(Lon);
Z=(N*(1-e2)+H).*sin(Lat);