围绕上一篇，同一椭球基准内不同坐标系下坐标转换的相关算法，此篇主要分享的是这些算法的实现
采用代码为C#，转换精度都做了验证，都在小数点后两位
大地球面坐标与大地空间直角坐标互转
【坐标正转】 经度、纬度、高程(LBH)转到XYZ
/// <summary>
/// 大地坐标转空间直角坐标
/// </summary>
/// <param name="geoCoordSystem"></param>
/// <param name="originalPoint">纬度、经度、大地高。</param>
/// <param name="reslutPoint"></param>
public void GeodeticToSpace(GeoCoordSystem geoCoordSystem,
GeoPoint originalPoint,out Point3d reslutPoint)
{
double B, L, H; //纬度、经度、大地高“十进制弧度”，大地高单位为:"米"
double a, b, f; //长半轴、短半轴、椭球反扁率
double e, ee, e2; //第一偏心率、第二偏心率，第一偏心率平方。
///参数示例
{
//a = 6378245.00;//椭球长半轴,单位米
//f = 1 / 298.30000000000001;//椭球扁率
//b = a * (1 - f);//椭球短半轴
//e = Math.Sqrt(a * a - b * b) / a;//椭球第一偏心率
//ee = Math.Sqrt(a * a - b * b) / b;//椭球第二偏心率
//e2 = e * e; //第一偏心率的平方。
}
B = UnitConvert.DegreesToRedians(originalPoint.N);
L = UnitConvert.DegreesToRedians(originalPoint.E);
H = originalPoint.H;
a = geoCoordSystem.SpheroidDatum.Semiminor_Axis;
e = geoCoordSystem.SpheroidDatum.Spheroid_E;
e2 = e * e;
double N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B)); //子午圈曲率半径
//double X = (N + H) * Math.Cos(B) * Math.Cos(L);
//double Y = (N + H) * Math.Cos(B) * Math.Sin(L);
//double Z = (N * (1 - e2) + H) * Math.Sin(B);
reslutPoint.X = (N + H) * Math.Cos(B) * Math.Cos(L);
reslutPoint.Y = (N + H) * Math.Cos(B) * Math.Sin(L);
reslutPoint.Z = (N * (1 - e2) + H) * Math.Sin(B);
}【坐标反转】 XYZ转经度、纬度、高程(LBH)，空间直角坐标转经纬度坐标，迭代法的实现
/// <summary>
/// 空间直角坐标转大地坐标
/// </summary>
/// <param name="geoCoordSystem">参考坐标系</param>
/// <param name="originalPoint">空间直角坐标</param>
/// <param name="reslutPoint">返回结果：大地坐标</param>
public void SpaceRectangularToGeodetic(GeoCoordSystem geoCoordSystem
, Point3d originalPoint, out GeoPoint reslutPoint)
{
double B, L, H; //经度、纬度单位此处为:“十进制弧度”，大地高单位为:"米"
double a, b, f; //长半轴、短半轴、椭球反扁率
double e, ee, e2; //第一偏心率、第二偏心率，第一偏心率平方。
double X, Y, Z; //空间直角坐标点的，X,Y,Z
double Rho= 206264.806247096363;
var Error = 0.00001;
///参数示例
{
//a = 6378245.00;//椭球长半轴,单位米
//f = 1 / 298.30000000000001;//椭球扁率
//b = a * (1 - f);//椭球短半轴
//e = Math.Sqrt(a * a - b * b) / a;//椭球第一偏心率
//ee = Math.Sqrt(a * a - b * b) / b;//椭球第二偏心率
//e2 = e * e; //第一偏心率的平方。
}
X = originalPoint.X;
Y = originalPoint.Y;
Z = originalPoint.Z;
a = geoCoordSystem.SpheroidDatum.Semiminor_Axis;
e = geoCoordSystem.SpheroidDatum.Spheroid_E;
e2 = e * e;
double B1 = Math.Atan(Z / Math.Sqrt(X * X + Y * Y));
double N1, B2;
while (true)
{
N1 = a / Math.Sqrt(1 - e2 * Math.Sin(B1) * Math.Sin(B1)); //子午圈曲率半径
B2 = Math.Atan((Z + N1 * e2 * Math.Sin(B1)) / Math.Sqrt(X * X + Y * Y));
if (Math.Abs(B2 - B1) * Rho <= Error)
break;
B1 = B2;
}
B = B2;
L = Math.Acos(X / Math.Sqrt(X * X + Y * Y));
if (Y < 0)
{
L = -L;
}
double N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B)); //子午圈曲率半径
H= Math.Sqrt(X * X + Y * Y) / Math.Cos(B) - N;
reslutPoint.N = UnitConvert.RediansToDegrees(B);
reslutPoint.E = UnitConvert.RediansToDegrees(L);
reslutPoint.H =H;
}大地球面坐标与大地投影坐标互转-以高斯投影坐标转换为例
【坐标正算】 经度、纬度(LB)转到x,y
/// <summary>
/// 大地坐标转高斯坐标
/// </summary>
/// <param name="pGaussKrugerCoordSystem"></param>
/// <param name="originalCooords"></param>
/// <param name="targetCooords"></param>
public void GeodeticToGauss(GaussKrugerCoordSystem pGaussKrugerCoordSystem,
GeoPoint originalCooords, out GeoPoint targetCooords)
{
double B, L, H; //纬度、经度、大地高。单位为十进制度。
double a, b, f; //长半轴、短半轴、椭球反扁率
double e, ee, e2; //第一偏心率、第二偏心率，第一偏心率平方。
double x, y, z; //高斯坐标X,Y,Z
double L0; //中央子午线（中央经线）
{
//double x, y;
//double B = 40.00; //纬度 单位度
//double L = 113.00; //经度 单位 度
//B = B * Math.PI / 180; ///0.0174533; //转换为弧度
//L = L * Math.PI / 180; //0.0174533; //转换为弧度
//double a = 6378245.00;//椭球长半轴,单位米
//double f = 1 / 298.30000000000001;//椭球扁率
//double b = a * (1 - f);//椭球短半轴
//double e = Math.Sqrt(a * a - b * b) / a;//椭球第一偏心率
//double ee = Math.Sqrt(a * a - b * b) / b;//椭球第二偏心率
//double L0 = 114.00;//中央子午线 单位为度。
//L0 = L0 * Math.PI / 180; //0.0174533; //转换为弧度
//double e2 = e * e; //第一偏心率的平方。
}
B = UnitConvert.DegreesToRedians(originalCooords.N);
L = UnitConvert.DegreesToRedians(originalCooords.E);
H = UnitConvert.DegreesToRedians(originalCooords.H);
L0 = UnitConvert.DegreesToRedians(pGaussKrugerCoordSystem.Central_Meridian);
a = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Semiminor_Axis;
e = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Spheroid_E;
ee = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Spheroid_EE;
e2 = e * e;
double Sin_B = Math.Sin(B);
double Cos_B = Math.Cos(B);
//子午圈曲率半径
double N = a / Math.Sqrt(1 - Math.Pow(e, 2) * Sin_B * Sin_B);
double ll = L - L0; //经差，到中央经线的差
double t2 = Math.Tan(B) * Math.Tan(B);
double n2 = ee * ee * Cos_B * Cos_B;
double pp = 1.0;// (180 / Math.PI) * 3600; //
double m0 = a * (1 - e * e);
double m2 = 3 * e2 * m0 / 2.0;
double m4 = 5 * e2 * m2 / 4.0;
double m6 = 7 * e2 * m4 / 6.0;
double m8 = 9 * e2 * m6 / 8.0;
double a0 = m0 + m2 / 2.0 + 3 * m4 / 8.0 + 5 * m6 / 16.0 + 35 * m8 / 128.0;
double a2 = m2 / 2.0 + m4 / 2.0 + 15 * m6 / 32.0 + 7 * m8 / 16.0;
double a4 = m4 / 8.0 + 3 * m6 / 16.0 + 7 * m8 / 32.0;
double a6 = m6 / 32.0 + m8 / 16.0;
double a8 = m8 / 128.0;
//中央子午线弧长。
double X = a0 * B - Sin_B * Cos_B * ((a2 - a4 + a6)
+ (2 * a4 - 16 * a6 / 3.0) * Sin_B * Sin_B
+ 16 * a6 * Math.Pow(Sin_B, 4) / 3.0);
double x1 = N * Sin_B * Cos_B * ll * ll / (2 * pp * pp);
double x2 = N * Sin_B * Cos_B * Cos_B * Cos_B
*(5 - t2 + 9 * n2 + 4 * n2 * n2) * Math.Pow(ll, 4)
/ (24 * Math.Pow(pp, 4));
double x3 = N * Sin_B * Math.Pow(Cos_B, 5)
* (61 - 58 * t2 + t2 * t2) * Math.Pow(ll, 6)
/ (720 * Math.Pow(pp, 6));
x = X + x1 + x2 + x3;
double y1 = N * Cos_B * ll / pp;
double y2 = N * Math.Pow(Cos_B, 3) * (1 - t2 + n2) * Math.Pow(ll, 3)
/ (6 * Math.Pow(pp, 3));
double y3 = N * Math.Pow(Cos_B, 5)
* (5 - 18 * t2 + t2 * t2 + 14 * n2 - 58 * n2 * t2)
* Math.Pow(ll, 5)
/ (120 * Math.Pow(pp, 5));
y = y1 + y2 + y3;
y = y + 500000;
targetCooords.N = x;
targetCooords.E = y;
targetCooords.H = 0;
}【坐标反算】 x,y转到经度、纬度(LB)
/// <summary>
/// 高斯坐标转大地坐标
/// </summary>
/// <param name="pGaussKrugerCoordSystem">高斯坐标参考</param>
/// <param name="originalCooords">高斯平面坐标点</param>
/// <param name="targetCooords">大地坐标点</param>
public void GaussToGeodetic(GaussKrugerCoordSystem pGaussKrugerCoordSystem,
GeoPoint originalCooord,out GeoPoint targetCooord)
{
double B, L, H; //纬度、经度、大地高。单位为十进制度。
double a, b, f; //长半轴、短半轴、椭球反扁率
double e, ee, e2; //第一偏心率、第二偏心率，第一偏心率平方。
double x, y, z; //高斯坐标X,Y,Z
double L0; //中央子午线（中央经线）
//double x = 4430486.4133464955, y = 414603.95129119331;
//y = y - 500000;
//double B = 0; //纬度 单位度, 返回结果
//double L = 0; //经度 单位 度，返回结果
//double a = 6378245.00; //椭球长半轴,单位米
//double f = 1 / 298.30000000000001; //椭球扁率
//double b = a * (1 - f); //椭球短半轴
//double e = Math.Sqrt(a * a - b * b) / a;//椭球第一偏心率
//double ee = Math.Sqrt(a * a - b * b) / b;//椭球第二偏心率
//double L0 = 114.00;//中央子午线 单位为度。
//L0 = L0 * Math.PI / 180; //0.0174533; //转换为弧度
//double e2 = e * e; //第一偏心率的平方。
x = originalCooord.N- pGaussKrugerCoordSystem.False_Northing;
y = originalCooord.E- pGaussKrugerCoordSystem.False_Easting;
z = originalCooord.H;
L0 = UnitConvert.DegreesToRedians(pGaussKrugerCoordSystem.Central_Meridian);
a = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Semiminor_Axis;
e = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Spheroid_E;
ee = pGaussKrugerCoordSystem.GeoCoordSystem.SpheroidDatum.Spheroid_EE;
e2 = e * e;
double m0 = a * (1 - e * e);
double m2 = 3 * e2 * m0 / 2.0;
double m4 = 5 * e2 * m2 / 4.0;
double m6 = 7 * e2 * m4 / 6.0;
double m8 = 9 * e2 * m6 / 8.0;
double a0 = m0 + m2 / 2.0 + 3 * m4 / 8.0 + 5 * m6 / 16.0 + 35 * m8 / 128.0;
double a2 = m2 / 2.0 + m4 / 2.0 + 15 * m6 / 32.0 + 7 * m8 / 16.0;
double a4 = m4 / 8.0 + 3 * m6 / 16.0 + 7 * m8 / 32.0;
double a6 = m6 / 32.0 + m8 / 16.0;
double a8 = m8 / 128.0;
double Bf = 0;//低点纬度
double Bfi = x / a0;
double Bfi1 = Bfi;
int m = 0;
do
{
Bfi = Bfi1;
Bfi1 = (x - (-a2 * Math.Sin(2 * Bfi) / 2
+ a4 * Math.Sin(4 * Bfi) / 4 - a6 * Math.Sin(6 * Bfi) / 6
+ a8 * Math.Sin(8 * Bfi) / 8))
/ a0;
m++;
} while (Math.Abs(Bfi1 - Bfi) > 1e-15);
Bf = Bfi;
double Nf = a * Math.Pow((1 - Math.Sin(Bf) * Math.Sin(Bf) * e * e), -0.5);
double Mf = a * (1 - e * e) * Math.Pow(1 - Math.Sin(Bf)
* Math.Sin(Bf) * e * e, -1.5);
double tf = Math.Tan(Bf);
double nf = ee * Math.Cos(Bf);
double W = 180 / Math.PI; //弧度转度。
B = Bf - tf * y * y / (2 * Mf * Nf) +
tf * (5 + 3 * tf * tf + nf * nf
- 9 * nf * nf * tf * tf) * Math.Pow(y, 4) / (24 * Mf * Math.Pow(Nf, 3))
-tf * (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * Math.Pow(y, 6)
/ (720 * Mf * Math.Pow(Nf, 5));
L = L0
+ y / Nf / Math.Cos(Bf) - Math.Pow(y, 3)
* (1 + 2 * tf * tf + nf * nf) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf))
+ Math.Pow(y, 5) * (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4)
+ 6 * nf * nf + 8 * nf * nf * tf * tf)
/ (120 * Math.Pow(Nf, 5) * Math.Cos(Bf));
targetCooord.N = B;
targetCooord.E = L;
B = B * 180 / Math.PI;
L = L * 180 / Math.PI;
targetCooord.N = B;
targetCooord.E = L;
targetCooord.H = 0;
}