序:网上很多算法,基本都是错的,还得去找论文才行。原理就不讲了,我会贴出论文引用。
1.转换过程
(1)确定站心直角坐标系的原点(X0,Y0,Z0)=(0,0,0)的经纬高(B0,L0,H0)作为基准点。
(2)将基准点(B0,L0,H0)转成地心坐标。
(3)将站心直角坐标系的目标点(X1,Y1,Z1)转成地心坐标(eX1,eY1,eZ1)。
(4)将地心坐标(eX1,eY1,eZ1)转成经纬高(B1,L1,H1)。
2.C#代码
/// <summary>
/// 使用WGS84椭球体
/// </summary>
public class CoorTransform
{
#region C#实现方法/三角函数及迭代法
/// <summary>
/// 地心坐标转经纬高
/// </summary>
/// <param name="X"></param>
/// <param name="Y"></param>
/// <param name="Z"></param>
/// <param name="aAxis"></param>
/// <param name="bAxis"></param>
/// <returns></returns>
public static Vector3D XYZtoBLH(double X, double Y, double Z, double aAxis = 6378137, double bAxis = 6356752.3142)
{
double targetH, targetB, targetL;
double e1 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(aAxis, 2);
double e2 = (Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(bAxis, 2);
double S = Math.Sqrt(Math.Pow(X, 2) + Math.Pow(Y, 2));
double cosL = X / S;
double B = 0;
double L = 0;
L = Math.Acos(cosL);
L = Math.Abs(L);
double tanB = Z / S;
B = Math.Atan(tanB);
double c = aAxis * aAxis / bAxis;
double preB0 = 0.0;
double ll = 0.0;
double N = 0.0;
//迭代计算纬度
do
{
preB0 = B;
ll = Math.Pow(Math.Cos(B), 2) * e2;
N = c / Math.Sqrt(1 + ll);
tanB = (Z + N * e1 * Math.Sin(B)) / S;
B = Math.Atan(tanB);
} while (Math.Abs(preB0 - B) >= 0.0000000001);
ll = Math.Pow(Math.Cos(B), 2) * e2;
N = c / Math.Sqrt(1 + ll);
targetH = S / Math.Cos(B) - N;// * (1 - e1);
targetB = B * 180 / Math.PI;
targetL = L * 180 / Math.PI;
return new Vector3D(targetB, targetL, targetH);
}
/// <summary>
/// 经纬高转地心坐标
/// </summary>
/// <param name="B"></param>
/// <param name="L"></param>
/// <param name="H"></param>
/// <param name="aAxis"></param>
/// <param name="bAxis"></param>
/// <returns></returns>
public static Vector3D BLHtoXYZ(double B, double L, double H, double aAxis = 6378137, double bAxis = 6356752.3142)
{
double targetX, targetY, targetZ;
double dblD2R = Math.PI / 180.0;
double f = 1 / 298.257223563;
double e2 = 2 * f - f * f;
//double e1 = Math.Sqrt(Math.Pow(aAxis, 2) - Math.Pow(bAxis, 2)) / Math.Pow(aAxis, 2);
double W = Math.Sqrt(1 - e2 * Math.Sin(B * dblD2R) * Math.Sin(B * dblD2R));
double N = aAxis / W;
//double N = aAxis / Math.Sqrt(1.0 - Math.Pow(e1, 2) * Math.Pow(Math.Sin(B * dblD2R), 2));
targetX = (N + H) * Math.Cos(B * dblD2R) * Math.Cos(L * dblD2R);
targetY = (N + H) * Math.Cos(B * dblD2R) * Math.Sin(L * dblD2R);
targetZ = (N * (1.0 - e2) + H) * Math.Sin(B * dblD2R);
return new Vector3D(targetX, targetY, targetZ);
}
/// <summary>
/// 计算站心直角坐标系中的相对坐标在地心坐标系中的坐标
/// </summary>
/// <param name="X0">站心原点的地心坐标X</param>
/// <param name="Y0">站心原点的地心坐标Y</param>
/// <param name="Z0">站心原点的地心坐标Z</param>
/// <param name="B0">站心原点的纬度</param>
/// <param name="L0">站心原点的经度</param>
/// <param name="H0">站心原点的高度</param>
/// <param name="rX">相对站心原点的目标坐标X</param>
/// <param name="rY">相对站心原点的目标坐标Y</param>
/// <param name="rZ">相对站心原点的目标坐标Z</param>
/// <returns></returns>
public static Vector3D CalRelaXYZBasedOnBLH(double X0,double Y0,double Z0,double B0,double L0,double H0,
double rX,double rY,double rZ)
{
double dblD2R = Math.PI / 180.0;
var AX = -Math.Sin(L0 * dblD2R) * rX + -Math.Sin(B0 * dblD2R) * Math.Cos(L0 * dblD2R) * rY + Math.Cos(B0 * dblD2R) * Math.Cos(L0 * dblD2R) * rZ;
var AY = -Math.Cos(L0 * dblD2R) * rX + -Math.Sin(B0 * dblD2R) * Math.Sin(L0 * dblD2R) * rY + Math.Cos(B0 * dblD2R) * Math.Sin(L0 * dblD2R) * rZ;
var AZ = 0 + Math.Cos(B0 * dblD2R) * rY + Math.Sin(B0 * dblD2R) * rZ;
var relaX = AX + X0;
var relaY = AY + Y0;
var relaZ = AZ + Z0;
return new Vector3D(relaX, relaY, relaZ);
}
#endregion
}
3.转换代码,为了直观变量名使用中文
double 基准点纬度 = 34, 基准点经度 = 112, 基准点高 = 300;
double 目标站心坐标X = 0, 目标站心坐标Y = 100000, 目标站心坐标Z = 8000;
var 站心原点的地心坐标 = CoorTransform.BLHtoXYZ(基准点纬度, 基准点经度, 基准点高);
var 目标的地心坐标 = CoorTransform.CalRelaXYZBasedOnBLH(站心原点的地心坐标.X, 站心原点的地心坐标.Y, 站心原点的地心坐标.Z, 基准点纬度, 基准点经度, 基准点高, 目标站心坐标X, 目标站心坐标Y, 目标站心坐标Z);
var 目标的纬经高 = CoorTransform.XYZtoBLH(目标的地心坐标.X, 目标的地心坐标.Y, 目标的地心坐标.Z);
4.结果
[注:]
1.直角坐标单位皆为米。
2.本文的站心直角坐标为东、北、天,即X+ => 东,Y+ => 北,Z+ => 天。如果使用北、天、东注意下转换。
引用
[1]凌震莹.大地坐标系与站心地平直角坐标系的坐标转换[J].声学与电子工程,2009(04):31-34.