大地坐标与空间直角坐标转换

一、椭球
1、常用椭球参数

	克拉索夫斯基椭球	1975年国际椭球	WGS-84椭球	CGCS2000(2000中国大地坐标系)
a	6378245(m)	6378140(m)	6378137(m)	6378137(m)
b	6356863.0187730473(m)	6356755.2881575287(m)	6356752.3142 (m)	6356752.3141(m)
c	6399698.9017827110(m)	6399596.6519880105(m)	6399593.6258(m)	6399593.6259(m)
α	1/298.3	1/298.257	1/298.257223563	1/298.257222101
e2	0.006693421622966	0.006694384999588	0.00669437999013	0.00669438002290
e’2	0.006738525414683	0.006739501819473	0.00673949674227	0.00673949677548

决定旋转椭球的大小和形状只要知道以上参数中的两个就可（至少包含一个长度元素）

我国的北京54坐标系是基于克拉索夫斯基椭球建立的，
西安80坐标系是基于1975年国际椭球建立，
GPS定位系统是基于WGS-84椭球，

CGCS2000坐标系是全球地心坐标系在我国的体现，坐标系原点是地球质心（包括大气和海洋），
Z轴指向国际时间局1984年定义的CTP，X轴为起始子午面与通过原点且正交于Z轴的赤道面，构成右手地心地固直角坐标系

2、
C#构建椭球类

class Ellipsoid
    {
        //长半轴
        public double a;
        //短半轴
        public double b;
        //扁率
        public double f;
        //第一偏心率平方
        public double ec;
        //第二偏心率平方
        public double ecc;
        /// <summary>
        /// 
        /// </summary>
        /// <param name="a为长半轴"></param>
        /// <param name="Invf为扁率"></param>
        public Ellipsoid(double a,double Invf)
        {
            this.a = a;
            f = 1 / Invf;
            evaluate();
        }
        public Ellipsoid()
        {
            a = 6378137.000;
            f = 1 / 298.3;
            evaluate();
        }
        void evaluate()
        {
            double b = a - a * f;
            ec = 1 - Math.Pow(b, 2) / Math.Pow(a, 2);
            ecc= ec / (1-ec);
        }
    }

二、大地坐标与空间直角坐标转换

1、转换类

class Transform
    {
        public double a;
        public double ec;
        public double ecc;
        public Transform(Ellipsoid ellipsoid)
        {
            a = ellipsoid.a;
            ec = ellipsoid.ec;
            ecc = ellipsoid.ecc;
        }
        public PointXYZ BLHToXYZ(PointBLH  bLH)
        {
            double B = bLH.B;
            double L = bLH.L;
            double H = bLH.H;
            //辅助计算公式
            double W = Math.Sqrt(1 - ec * Math.Pow(Math.Sin(B), 2));
            double N = a / W;
            PointXYZ xYZ = new PointXYZ();
            xYZ.X = (N + H) * Math.Cos(B) * Math.Cos(L);
            xYZ.Y = (N + H) * Math.Cos(B) * Math.Sin(L);
            xYZ.Z = (N * (1 - ec) + H) * Math.Sin(B);
            return xYZ;
        }
        
        public PointBLH XYZToBLH(PointXYZ xYZ)
        {
            double X = xYZ.X;
            double Y = xYZ.Y;
            double Z = xYZ.Z;
            PointBLH bLH = new PointBLH();
            bLH.L = Math.Atan(Y / X);
            if (X < 0)
                bLH.L += Math.PI;
            var r = Math.Sqrt(X * X + Y * Y);
            var B1 = Math.Atan(Z / r);
            double B2;
            while (true)
            {
                var W1 = Math.Sqrt(1 - ec * (Math.Sin(B1) * Math.Sin(B1)));
                var N1 = a / W1;
                B2 = Math.Atan((Z + N1 * ec * Math.Sin(B1)) / r);
                if (Math.Abs(B2 - B1) <= 0.0000000001)
                    break;
                B1 = B2;
            }
            bLH.B = B2;
            var W = Math.Sqrt(1 - ec * (Math.Sin(B2) * Math.Sin(B2)));
            var N = a / W;
            bLH.H = r / Math.Cos(B2) - N;
            return bLH;
        }
    }

把输入的点构造成Mypoint类

class PointBLH
    {
        public double B;
        public double L;
        public double H;
    }
    class PointXYZ
    {
        public double X;
        public double Y;
        public double Z;
    }

2、大地坐标系与空间之间坐标系转换公式

大地坐标换算到空间直角坐标
P点投影到椭球面上，大地高为H

空间直角坐标换算到大地坐标

解算大地纬度的时候需要进行迭代计算，
保证前后两次迭代<e-10时，可保证计算精度达到0.0001″